conformalVoronoiMeshCalcDualMesh.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2012-2016 OpenFOAM Foundation
    Copyright (C) 2018-2022 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 "conformalVoronoiMesh.H"
#include "motionSmoother.H"
#include "backgroundMeshDecomposition.H"
#include "polyMeshGeometry.H"
#include "indexedCellChecks.H"
#include "OBJstream.H"
#include "indexedCellOps.H"
#include "ListOps.H"
#include "DelaunayMeshTools.H"

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

void Foam::conformalVoronoiMesh::calcDualMesh
(
    pointField& points,
    labelList& boundaryPts,
    faceList& faces,
    labelList& owner,
    labelList& neighbour,
    wordList& patchNames,
    PtrList<dictionary>& patchDicts,
    pointField& cellCentres,
    labelList& cellToDelaunayVertex,
    labelListList& patchToDelaunayVertex,
    bitSet& boundaryFacesToRemove
)
{
    timeCheck("Start calcDualMesh");

    setVertexSizeAndAlignment();

    timeCheck("After setVertexSizeAndAlignment");

    indexDualVertices(points, boundaryPts);

    {
        Info<< nl << "Merging identical points" << endl;

        // There is no guarantee that a merge of close points is no-risk
        mergeIdenticalDualVertices(points, boundaryPts);
    }

    // Final dual face and owner neighbour construction

    timeCheck("Before createFacesOwnerNeighbourAndPatches");

    createFacesOwnerNeighbourAndPatches
    (
        points,
        faces,
        owner,
        neighbour,
        patchNames,
        patchDicts,
        patchToDelaunayVertex,  // from patch face to Delaunay vertex (slavePp)
        boundaryFacesToRemove,
        false
    );

    // deferredCollapseFaceSet(owner, neighbour, deferredCollapseFaces);

    cellCentres = DelaunayMeshTools::allPoints(*this);

    cellToDelaunayVertex = removeUnusedCells(owner, neighbour);

    cellCentres = pointField(cellCentres, cellToDelaunayVertex);

    removeUnusedPoints(faces, points, boundaryPts);

    timeCheck("End of calcDualMesh");
}


void Foam::conformalVoronoiMesh::calcTetMesh
(
    pointField& points,
    labelList& pointToDelaunayVertex,
    faceList& faces,
    labelList& owner,
    labelList& neighbour,
    wordList& patchNames,
    PtrList<dictionary>& patchDicts
)
{
    labelList vertexMap(number_of_vertices());

    label vertI = 0;

    points.setSize(number_of_vertices());
    pointToDelaunayVertex.setSize(number_of_vertices());

    for
    (
        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
        vit != finite_vertices_end();
        ++vit
    )
    {
        if (vit->internalPoint() || vit->boundaryPoint())
        {
            vertexMap[vit->index()] = vertI;
            points[vertI] = topoint(vit->point());
            pointToDelaunayVertex[vertI] = vit->index();
            vertI++;
        }
    }

    points.setSize(vertI);
    pointToDelaunayVertex.setSize(vertI);

    label celli = 0;

    for
    (
        Delaunay::Finite_cells_iterator cit = finite_cells_begin();
        cit != finite_cells_end();
        ++cit
    )
    {
        if (cit->internalOrBoundaryDualVertex())
        {
             cit->cellIndex() = celli++;
        }
        else
        {
            cit->cellIndex() = Cb::ctFar;
        }
    }

    patchNames = geometryToConformTo_.patchNames();

    patchNames.setSize(patchNames.size() + 1);

    patchNames[patchNames.size() - 1] = "foamyHexMesh_defaultPatch";

    label nPatches = patchNames.size();

    List<DynamicList<face>> patchFaces(nPatches);
    List<DynamicList<label>> patchOwners(nPatches);

    faces.setSize(number_of_finite_facets());

    owner.setSize(number_of_finite_facets());

    neighbour.setSize(number_of_finite_facets());

    label facei = 0;

    labelList verticesOnTriFace(3, label(-1));

    face newFace(verticesOnTriFace);

    for
    (
        Delaunay::Finite_facets_iterator fit = finite_facets_begin();
        fit != finite_facets_end();
        ++fit
    )
    {
        const Cell_handle c1(fit->first);
        const label oppositeVertex = fit->second;
        const Cell_handle c2(c1->neighbor(oppositeVertex));

        if (c1->hasFarPoint() && c2->hasFarPoint())
        {
            // Both tets are outside, skip
            continue;
        }

        label c1I = c1->cellIndex();
        label c2I = c2->cellIndex();

        label ownerCell = -1;
        label neighbourCell = -1;

        for (label i = 0; i < 3; i++)
        {
            verticesOnTriFace[i] = vertexMap
            [
                c1->vertex(vertex_triple_index(oppositeVertex, i))->index()
            ];
        }

        newFace = face(verticesOnTriFace);

        if (c1->hasFarPoint() || c2->hasFarPoint())
        {
            // Boundary face...
            if (c1->hasFarPoint())
            {
                //... with c1 outside
                ownerCell = c2I;
            }
            else
            {
                // ... with c2 outside
                ownerCell = c1I;

                reverse(newFace);
            }

            label patchIndex = geometryToConformTo_.findPatch
            (
                newFace.centre(points)
            );

            if (patchIndex == -1)
            {
                patchIndex = patchNames.size() - 1;

                WarningInFunction
                    << newFace.centre(points) << nl
                    << "did not find a surface patch. Adding to "
                    << patchNames[patchIndex]
                    << endl;
            }

            patchFaces[patchIndex].append(newFace);
            patchOwners[patchIndex].append(ownerCell);
        }
        else
        {
            // Internal face...
            if (c1I < c2I)
            {
                // ...with c1 as the ownerCell
                ownerCell = c1I;
                neighbourCell = c2I;

                reverse(newFace);
            }
            else
            {
                // ...with c2 as the ownerCell
                ownerCell = c2I;
                neighbourCell = c1I;
            }

            faces[facei] = newFace;
            owner[facei] = ownerCell;
            neighbour[facei] = neighbourCell;
            facei++;
        }
    }

    label nInternalFaces = facei;

    faces.setSize(nInternalFaces);
    owner.setSize(nInternalFaces);
    neighbour.setSize(nInternalFaces);

    sortFaces(faces, owner, neighbour);

//    bitSet boundaryFacesToRemove;
//    List<DynamicList<bool>> indirectPatchFace;
//
//    addPatches
//    (
//        nInternalFaces,
//        faces,
//        owner,
//        patchDicts,
//        boundaryFacesToRemove,
//        patchFaces,
//        patchOwners,
//        indirectPatchFace
//    );
}


void Foam::conformalVoronoiMesh::mergeIdenticalDualVertices
(
    const pointField& pts,
    labelList& boundaryPts
)
{
    // Assess close points to be merged

    label nPtsMerged = 0;
    label nPtsMergedSum = 0;

    do
    {
        Map<label> dualPtIndexMap;

        nPtsMerged = mergeIdenticalDualVertices
        (
            pts,
            dualPtIndexMap
        );

        reindexDualVertices(dualPtIndexMap, boundaryPts);

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

        nPtsMergedSum += nPtsMerged;

    } while (nPtsMerged > 0);

    if (nPtsMergedSum > 0)
    {
        Info<< "    Merged " << nPtsMergedSum << " points " << endl;
    }
}


Foam::label Foam::conformalVoronoiMesh::mergeIdenticalDualVertices
(
    const pointField& pts,
    Map<label>& dualPtIndexMap
) const
{
    label nPtsMerged = 0;

    for
    (
        Delaunay::Finite_facets_iterator fit = finite_facets_begin();
        fit != finite_facets_end();
        ++fit
    )
    {
        const Cell_handle c1(fit->first);
        const label oppositeVertex = fit->second;
        const Cell_handle c2(c1->neighbor(oppositeVertex));

        if (is_infinite(c1) || is_infinite(c2))
        {
            continue;
        }

        label& c1I = c1->cellIndex();
        label& c2I = c2->cellIndex();

        if ((c1I != c2I) && !c1->hasFarPoint() && !c2->hasFarPoint())
        {
            const Foam::point& p1 = pts[c1I];
            const Foam::point& p2 = pts[c2I];

            if (p1 == p2)
            {
//                if (c1->parallelDualVertex() || c2->parallelDualVertex())
//                {
//                    if (c1->vertexLowestProc() < c2->vertexLowestProc())
//                    {
//                        dualPtIndexMap.insert(c1I, c1I);
//                        dualPtIndexMap.insert(c2I, c1I);
//                    }
//                    else
//                    {
//                        dualPtIndexMap.insert(c1I, c2I);
//                        dualPtIndexMap.insert(c2I, c2I);
//                    }
//                }
                if (c1I < c2I)
                {
                    dualPtIndexMap.insert(c1I, c1I);
                    dualPtIndexMap.insert(c2I, c1I);
                }
                else
                {
                    dualPtIndexMap.insert(c1I, c2I);
                    dualPtIndexMap.insert(c2I, c2I);
                }

                nPtsMerged++;
            }
        }
    }

    if (debug)
    {
        Info<< "mergeIdenticalDualVertices:" << endl
            << "    zero-length edges     : "
            << returnReduce(nPtsMerged, sumOp<label>()) << endl
            << endl;
    }

    return nPtsMerged;
}


//void Foam::conformalVoronoiMesh::smoothSurface
//(
//    pointField& pts,
//    const labelList& boundaryPts
//)
//{
//    label nCollapsedFaces = 0;
//
//    label iterI = 0;
//
//    do
//    {
//        Map<label> dualPtIndexMap;
//
//        nCollapsedFaces = smoothSurfaceDualFaces
//        (
//            pts,
//            boundaryPts,
//            dualPtIndexMap
//        );
//
//        reduce(nCollapsedFaces, sumOp<label>());
//
//        reindexDualVertices(dualPtIndexMap);
//
//        mergeIdenticalDualVertices(pts, boundaryPts);
//
//        if (nCollapsedFaces > 0)
//        {
//            Info<< "    Collapsed " << nCollapsedFaces << " boundary faces"
//                << endl;
//        }
//
//        if (++iterI > foamyHexMeshControls().maxCollapseIterations())
//        {
//            Info<< "    maxCollapseIterations reached, stopping collapse"
//                << endl;
//
//            break;
//        }
//
//    } while (nCollapsedFaces > 0);
//
//    // Force all points of boundary faces to be on the surface
//}
//
//
//Foam::label Foam::conformalVoronoiMesh::smoothSurfaceDualFaces
//(
//    pointField& pts,
//    const labelList& boundaryPts,
//    Map<label>& dualPtIndexMap
//) const
//{
//    label nCollapsedFaces = 0;
//
//    const scalar cosPerpendicularToleranceAngle = cos
//    (
//        degToRad(foamyHexMeshControls().surfaceStepFaceAngle())
//    );
//
//    for
//    (
//        Delaunay::Finite_edges_iterator eit = finite_edges_begin();
//        eit != finite_edges_end();
//        ++eit
//    )
//    {
//        Cell_circulator ccStart = incident_cells(*eit);
//        Cell_circulator cc = ccStart;
//
//        bool skipFace = false;
//
//        do
//        {
//            if (dualPtIndexMap.found(cc->cellIndex()))
//            {
//                // One of the points of this face has already been
//                // collapsed this sweep, leave for next sweep
//
//                skipFace = true;
//
//                break;
//            }
//
//        } while (++cc != ccStart);
//
//        if (skipFace)
//        {
//            continue;
//        }
//
//        if (isBoundaryDualFace(eit))
//        {
//            face dualFace = buildDualFace(eit);
//
//            if (dualFace.size() < 3)
//            {
//                // This face has been collapsed already
//                continue;
//            }
//
//            label maxFC = maxFilterCount(eit);
//
//            if (maxFC > foamyHexMeshControls().filterCountSkipThreshold())
//            {
//                // A vertex on this face has been limited too many
//                // times, skip
//                continue;
//            }
//
//            Cell_handle c = eit->first;
//            Vertex_handle vA = c->vertex(eit->second);
//            Vertex_handle vB = c->vertex(eit->third);
//
//            if
//            (
//                vA->internalBoundaryPoint() && vA->surfacePoint()
//             && vB->externalBoundaryPoint() && vB->surfacePoint()
//            )
//            {
//                if (vA->index() == vB->index() - 1)
//                {
//                    continue;
//                }
//            }
//            else if
//            (
//                vA->externalBoundaryPoint() && vA->surfacePoint()
//             && vB->internalBoundaryPoint() && vB->surfacePoint()
//            )
//            {
//                if (vA->index() == vB->index() + 1)
//                {
//                    continue;
//                }
//            }
//
//
//            if ((faceNormal & surfaceNormal) < cosPerpendicularToleranceAngle)
//            {
//                scalar targetFaceSize = averageAnyCellSize(vA, vB);
//
//                // Selecting faces to collapse based on angle to
//                // surface, so set collapseSizeLimitCoeff to GREAT to
//                // allow collapse of all faces
//
//                faceCollapseMode mode = collapseFace
//                (
//                    dualFace,
//                    pts,
//                    boundaryPts,
//                    dualPtIndexMap,
//                    targetFaceSize,
//                    GREAT,
//                    maxFC
//                );
//
//                if (mode == fcmPoint || mode == fcmEdge)
//                {
//                    nCollapsedFaces++;
//                }
//            }
//        }
//    }
//
//    return nCollapsedFaces;
//}


void Foam::conformalVoronoiMesh::deferredCollapseFaceSet
(
    labelList& owner,
    labelList& neighbour,
    const labelPairHashSet& deferredCollapseFaces
) const
{
    DynamicList<label> faceLabels;

    forAll(neighbour, nI)
    {
        if (deferredCollapseFaces.found(Pair<label>(owner[nI], neighbour[nI])))
        {
            faceLabels.append(nI);
        }
    }

    Pout<< "facesToCollapse" << nl << faceLabels << endl;
}


Foam::autoPtr<Foam::polyMesh>
Foam::conformalVoronoiMesh::createPolyMeshFromPoints
(
    const pointField& pts
) const
{
    faceList faces;
    labelList owner;
    labelList neighbour;
    wordList patchNames;
    PtrList<dictionary> patchDicts;
    pointField cellCentres;
    labelListList patchToDelaunayVertex;
    bitSet boundaryFacesToRemove;

    timeCheck("Start of checkPolyMeshQuality");

    Info<< nl << "Creating polyMesh to assess quality" << endl;

    createFacesOwnerNeighbourAndPatches
    (
        pts,
        faces,
        owner,
        neighbour,
        patchNames,
        patchDicts,
        patchToDelaunayVertex,
        boundaryFacesToRemove,
        false
    );

    cellCentres = DelaunayMeshTools::allPoints(*this);

    labelList cellToDelaunayVertex(removeUnusedCells(owner, neighbour));
    cellCentres = pointField(cellCentres, cellToDelaunayVertex);

    auto meshPtr = autoPtr<polyMesh>::New
    (
        IOobject
        (
            "foamyHexMesh_temporary",
            runTime_.timeName(),
            runTime_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        pointField(pts),  // Copy of points
        std::move(faces),
        std::move(owner),
        std::move(neighbour)
    );

    polyMesh& pMesh = meshPtr();

    List<polyPatch*> patches(patchNames.size());

    label nValidPatches = 0;

    forAll(patches, p)
    {
        label nPatchFaces = patchDicts[p].get<label>("nFaces");

        if
        (
            patchDicts.set(p)
         && patchDicts[p].get<word>("type") == processorPolyPatch::typeName
        )
        {
            // Do not create empty processor patches
            if (nPatchFaces)
            {
                patchDicts[p].set("transform", "coincidentFullMatch");

                patches[nValidPatches] = new processorPolyPatch
                (
                    patchNames[p],
                    patchDicts[p],
                    nValidPatches,
                    pMesh.boundaryMesh(),
                    processorPolyPatch::typeName
                );

                nValidPatches++;
            }
        }
        else
        {
            // Check that the patch is not empty on every processor

            if (returnReduceOr(nPatchFaces))
            {
                patches[nValidPatches] = polyPatch::New
                (
                    patchNames[p],
                    patchDicts[p],
                    nValidPatches,
                    pMesh.boundaryMesh()
                ).ptr();

                nValidPatches++;
            }
        }
    }

    patches.setSize(nValidPatches);

    pMesh.addPatches(patches);

    return meshPtr;
}


void Foam::conformalVoronoiMesh::checkCellSizing()
{
    Info<< "Checking cell sizes..."<< endl;

    timeCheck("Start of Cell Sizing");

    labelList boundaryPts(number_of_finite_cells(), internal);
    pointField ptsField;

    indexDualVertices(ptsField, boundaryPts);

    // Merge close dual vertices.
    mergeIdenticalDualVertices(ptsField, boundaryPts);

    autoPtr<polyMesh> meshPtr = createPolyMeshFromPoints(ptsField);
    const polyMesh& pMesh = meshPtr();

    //pMesh.write();

    // Find cells with poor quality
    DynamicList<label> checkFaces(identity(pMesh.nFaces()));
    labelHashSet wrongFaces(pMesh.nFaces()/100);

    Info<< "Running checkMesh on mesh with " << pMesh.nCells()
        << " cells "<< endl;

    const dictionary& dict
        = foamyHexMeshControls().foamyHexMeshDict();

    const dictionary& meshQualityDict
        = dict.subDict("meshQualityControls");

    const scalar maxNonOrtho =
        meshQualityDict.get<scalar>("maxNonOrtho", keyType::REGEX_RECURSIVE);

    label nWrongFaces = 0;

    if (maxNonOrtho < 180.0 - SMALL)
    {
        polyMeshGeometry::checkFaceDotProduct
        (
            false,
            maxNonOrtho,
            pMesh,
            pMesh.cellCentres(),
            pMesh.faceAreas(),
            checkFaces,
            List<labelPair>(),
            &wrongFaces
        );

        label nNonOrthogonal = returnReduce(wrongFaces.size(), sumOp<label>());

        Info<< "    non-orthogonality > " << maxNonOrtho
            << " degrees : " << nNonOrthogonal << endl;

        nWrongFaces += nNonOrthogonal;
    }

    labelHashSet protrudingCells = findOffsetPatchFaces(pMesh, 0.25);

    label nProtrudingCells = protrudingCells.size();

    Info<< "    protruding/intruding cells : " << nProtrudingCells << endl;

    nWrongFaces += nProtrudingCells;

//    motionSmoother::checkMesh
//    (
//        false,
//        pMesh,
//        meshQualityDict,
//        checkFaces,
//        wrongFaces
//    );

    Info<< "    Found total of " << nWrongFaces << " bad faces" << endl;

    {
        labelHashSet cellsToResizeMap(pMesh.nFaces()/100);

        // Find cells that are attached to the faces in wrongFaces.
        for (const label facei : wrongFaces)
        {
            const label faceOwner = pMesh.faceOwner()[facei];
            const label faceNeighbour = pMesh.faceNeighbour()[facei];

            cellsToResizeMap.insert(faceOwner);
            cellsToResizeMap.insert(faceNeighbour);
        }

        cellsToResizeMap += protrudingCells;

        pointField cellsToResize(cellsToResizeMap.size());

        label count = 0;
        for (label celli = 0; celli < pMesh.nCells(); ++celli)
        {
            if (cellsToResizeMap.found(celli))
            {
                cellsToResize[count++] = pMesh.cellCentres()[celli];
            }
        }

        Info<< "    DISABLED: Automatically re-sizing " << cellsToResize.size()
            << " cells that are attached to the bad faces: " << endl;

        //cellSizeControl_.setCellSizes(cellsToResize);
    }

    timeCheck("End of Cell Sizing");

    Info<< "Finished checking cell sizes"<< endl;
}


Foam::labelHashSet Foam::conformalVoronoiMesh::findOffsetPatchFaces
(
    const polyMesh& mesh,
    const scalar allowedOffset
) const
{
    timeCheck("Start findRemainingProtrusionSet");

    const polyBoundaryMesh& patches = mesh.boundaryMesh();

    cellSet offsetBoundaryCells
    (
        mesh,
        "foamyHexMesh_protrudingCells",
        mesh.nCells()/1000
    );

    forAll(patches, patchi)
    {
        const polyPatch& patch = patches[patchi];

        const faceList& localFaces = patch.localFaces();
        const pointField& localPoints = patch.localPoints();

        const labelList& fCell = patch.faceCells();

        forAll(localFaces, pLFI)
        {
            const face& f = localFaces[pLFI];

            const Foam::point& faceCentre = f.centre(localPoints);

            const scalar targetSize = targetCellSize(faceCentre);

            pointIndexHit pHit;
            label surfHit = -1;

            geometryToConformTo_.findSurfaceNearest
            (
                faceCentre,
                sqr(targetSize),
                pHit,
                surfHit
            );

            if
            (
                pHit.hit()
             && (pHit.point().dist(faceCentre) > allowedOffset*targetSize)
            )
            {
                offsetBoundaryCells.insert(fCell[pLFI]);
            }
        }
    }

    if (foamyHexMeshControls().objOutput())
    {
        offsetBoundaryCells.write();
    }

    return std::move(offsetBoundaryCells);
}


Foam::labelHashSet Foam::conformalVoronoiMesh::checkPolyMeshQuality
(
    const pointField& pts
) const
{
    autoPtr<polyMesh> meshPtr = createPolyMeshFromPoints(pts);
    polyMesh& pMesh = meshPtr();

    timeCheck("polyMesh created, checking quality");

    labelHashSet wrongFaces(pMesh.nFaces()/100);

    DynamicList<label> checkFaces(pMesh.nFaces());

    const vectorField& fAreas = pMesh.faceAreas();

    scalar faceAreaLimit = SMALL;

    forAll(fAreas, fI)
    {
        if (mag(fAreas[fI]) > faceAreaLimit)
        {
            checkFaces.append(fI);
        }
    }

    Info<< nl << "Excluding "
        << returnReduce(fAreas.size() - checkFaces.size(), sumOp<label>())
        << " faces from check, < " << faceAreaLimit << " area" << endl;

    const dictionary& dict
        = foamyHexMeshControls().foamyHexMeshDict();

    const dictionary& meshQualityDict
        = dict.subDict("meshQualityControls");

    motionSmoother::checkMesh
    (
        false,
        pMesh,
        meshQualityDict,
        checkFaces,
        wrongFaces
    );

    {
        // Check for cells with more than 1 but fewer than 4 faces
        label nInvalidPolyhedra = 0;

        const cellList& cells = pMesh.cells();

        forAll(cells, cI)
        {
            if (cells[cI].size() < 4 && cells[cI].size() > 0)
            {
                // Pout<< "cell " << cI << " " << cells[cI]
                //     << " has " << cells[cI].size() << " faces."
                //     << endl;

                nInvalidPolyhedra++;

                wrongFaces.insert(cells[cI]);
            }
        }

        Info<< "    cells with more than 1 but fewer than 4 faces          : "
            << returnReduce(nInvalidPolyhedra, sumOp<label>())
            << endl;

        // Check for cells with one internal face only

        labelList nInternalFaces(pMesh.nCells(), Zero);

        for (label fI = 0; fI < pMesh.nInternalFaces(); fI++)
        {
            nInternalFaces[pMesh.faceOwner()[fI]]++;
            nInternalFaces[pMesh.faceNeighbour()[fI]]++;
        }

        const polyBoundaryMesh& patches = pMesh.boundaryMesh();

        forAll(patches, patchi)
        {
            if (patches[patchi].coupled())
            {
                const labelUList& owners = patches[patchi].faceCells();

                forAll(owners, i)
                {
                    nInternalFaces[owners[i]]++;
                }
            }
        }

        label oneInternalFaceCells = 0;

        forAll(nInternalFaces, cI)
        {
            if (nInternalFaces[cI] <= 1)
            {
                oneInternalFaceCells++;
                wrongFaces.insert(cells[cI]);
            }
        }

        Info<< "    cells with with zero or one non-boundary face          : "
            << returnReduce(oneInternalFaceCells, sumOp<label>())
            << endl;
    }


    bitSet ptToBeLimited(pts.size(), false);

    for (const label facei : wrongFaces)
    {
        const face f = pMesh.faces()[facei];

        ptToBeLimited.set(f);
    }

    // // Limit connected cells

    // labelHashSet limitCells(pMesh.nCells()/100);

    // const labelListList& ptCells = pMesh.pointCells();

    // for (const label facei : wrongFaces)
    // {
    //     const face f = pMesh.faces()[facei];

    //     forAll(f, fPtI)
    //     {
    //         label ptI = f[fPtI];
    //         const labelList& pC = ptCells[ptI];
    //         limitCells.insert(pC);
    //     }
    // }

    // const labelListList& cellPts = pMesh.cellPoints();

    // for (const label celli : limitCells)
    // {
    //     const labelList& cP = cellPts[celli];

    //     ptToBeLimited.set(cP);
    // }


    // Apply Delaunay cell filterCounts and determine the maximum
    // overall filterCount

    label maxFilterCount = 0;

    for
    (
        Delaunay::Finite_cells_iterator cit = finite_cells_begin();
        cit != finite_cells_end();
        ++cit
    )
    {
        label cI = cit->cellIndex();

        if (cI >= 0)
        {
            if (ptToBeLimited.test(cI))
            {
                cit->filterCount()++;
            }

            if (cit->filterCount() > maxFilterCount)
            {
                maxFilterCount = cit->filterCount();
            }
        }
    }

    Info<< nl << "Maximum number of filter limits applied: "
        << returnReduce(maxFilterCount, maxOp<label>()) << endl;

    return wrongFaces;
}


Foam::label Foam::conformalVoronoiMesh::classifyBoundaryPoint
(
    Cell_handle cit
) const
{
    if (cit->boundaryDualVertex())
    {
        if (cit->featurePointDualVertex())
        {
            return featurePoint;
        }
        else if (cit->featureEdgeDualVertex())
        {
            return featureEdge;
        }
        else
        {
            return surface;
        }
    }
    else if (cit->baffleSurfaceDualVertex())
    {
        return surface;
    }
    else if (cit->baffleEdgeDualVertex())
    {
        return featureEdge;
    }
    else
    {
        return internal;
    }
}


void Foam::conformalVoronoiMesh::indexDualVertices
(
    pointField& pts,
    labelList& boundaryPts
)
{
    // Indexing Delaunay cells, which are the dual vertices

    this->resetCellCount();

    label nConstrainedVertices = 0;
    if (foamyHexMeshControls().guardFeaturePoints())
    {
        for
        (
            Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
            vit != finite_vertices_end();
            ++vit
        )
        {
            if (vit->constrained())
            {
                vit->index() = number_of_finite_cells() + nConstrainedVertices;
                nConstrainedVertices++;
            }
        }
    }

    pts.setSize(number_of_finite_cells() + nConstrainedVertices);
    boundaryPts.setSize
    (
        number_of_finite_cells() + nConstrainedVertices,
        internal
    );

    if (foamyHexMeshControls().guardFeaturePoints())
    {
        nConstrainedVertices = 0;
        for
        (
            Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
            vit != finite_vertices_end();
            ++vit
        )
        {
            if (vit->constrained())
            {
                pts[number_of_finite_cells() + nConstrainedVertices] =
                    topoint(vit->point());

                boundaryPts[number_of_finite_cells() + nConstrainedVertices] =
                    constrained;

                nConstrainedVertices++;
            }
        }
    }

    //OBJstream snapping1("snapToSurface1.obj");
    //OBJstream snapping2("snapToSurface2.obj");
    //OFstream tetToSnapTo("tetsToSnapTo.obj");

    for
    (
        Delaunay::Finite_cells_iterator cit = finite_cells_begin();
        cit != finite_cells_end();
        ++cit
    )
    {
//        if (tetrahedron(cit).volume() == 0)
//        {
//            Pout<< "ZERO VOLUME TET" << endl;
//            Pout<< cit->info();
//            Pout<< "Dual = " << cit->dual();
//        }

        if (!cit->hasFarPoint())
        {
            cit->cellIndex() = getNewCellIndex();

            // For nearly coplanar Delaunay cells that are present on different
            // processors the result of the circumcentre calculation depends on
            // the ordering of the vertices, so synchronise it across processors

            if (Pstream::parRun() && cit->parallelDualVertex())
            {
                typedef CGAL::Exact_predicates_exact_constructions_kernel Exact;
                typedef CGAL::Point_3<Exact> ExactPoint;

                List<labelPair> cellVerticesPair(4);
                List<ExactPoint> cellVertices(4);

                for (label vI = 0; vI < 4; ++vI)
                {
                    cellVerticesPair[vI] = labelPair
                    (
                        cit->vertex(vI)->procIndex(),
                        cit->vertex(vI)->index()
                    );

                    cellVertices[vI] = ExactPoint
                    (
                        cit->vertex(vI)->point().x(),
                        cit->vertex(vI)->point().y(),
                        cit->vertex(vI)->point().z()
                    );
                }

                // Sort the vertices so that they will be in the same order on
                // each processor
                labelList oldToNew(sortedOrder(cellVerticesPair));
                oldToNew = invert(oldToNew.size(), oldToNew);
                inplaceReorder(oldToNew, cellVertices);

                ExactPoint synchronisedDual = CGAL::circumcenter
                (
                    cellVertices[0],
                    cellVertices[1],
                    cellVertices[2],
                    cellVertices[3]
                );

                pts[cit->cellIndex()] = Foam::point
                (
                    CGAL::to_double(synchronisedDual.x()),
                    CGAL::to_double(synchronisedDual.y()),
                    CGAL::to_double(synchronisedDual.z())
                );
            }
            else
            {
                pts[cit->cellIndex()] = cit->dual();
            }

            // Feature point snapping
            if (foamyHexMeshControls().snapFeaturePoints())
            {
                if (cit->featurePointDualVertex())
                {
                    pointFromPoint dual = cit->dual();

                    pointIndexHit fpHit;
                    label featureHit;

                    // Find nearest feature point and compare
                    geometryToConformTo_.findFeaturePointNearest
                    (
                        dual,
                        sqr(targetCellSize(dual)),
                        fpHit,
                        featureHit
                    );

                    if (fpHit.hit())
                    {
                        if (debug)
                        {
                            Info<< "Dual        = " << dual << nl
                                << "    Nearest = " << fpHit.point() << endl;
                        }

                        pts[cit->cellIndex()] = fpHit.point();
                    }
                }
            }

//            {
//                // Snapping points far outside
//                if (cit->boundaryDualVertex() && !cit->parallelDualVertex())
//                {
//                    pointFromPoint dual = cit->dual();
//
//                    pointIndexHit hitInfo;
//                    label surfHit;
//
//                    // Find nearest surface point
//                    geometryToConformTo_.findSurfaceNearest
//                    (
//                        dual,
//                        sqr(targetCellSize(dual)),
//                        hitInfo,
//                        surfHit
//                    );
//
//                    if (!hitInfo.hit())
//                    {
//                        // Project dual to nearest point on tet
//
//                        tetPointRef tet
//                        (
//                            topoint(cit->vertex(0)->point()),
//                            topoint(cit->vertex(1)->point()),
//                            topoint(cit->vertex(2)->point()),
//                            topoint(cit->vertex(3)->point())
//                        );
//
//                        pointFromPoint nearestPointOnTet =
//                            tet.nearestPoint(dual).point();
//
//                        // Get nearest point on surface from tet.
//                        geometryToConformTo_.findSurfaceNearest
//                        (
//                            nearestPointOnTet,
//                            sqr(targetCellSize(nearestPointOnTet)),
//                            hitInfo,
//                            surfHit
//                        );
//
//                        vector snapDir = nearestPointOnTet - dual;
//                        snapDir /= mag(snapDir) + SMALL;
//
//                        drawDelaunayCell(tetToSnapTo, cit, offset);
//                        offset += 1;
//
//                        vectorField norm(1);
//                        allGeometry_[surfHit].getNormal
//                        (
//                            List<pointIndexHit>(1, hitInfo),
//                            norm
//                        );
//                        norm[0] /= mag(norm[0]) + SMALL;
//
//                        if
//                        (
//                            hitInfo.hit()
//                         && (mag(snapDir & norm[0]) > 0.5)
//                        )
//                        {
//                            snapping1.writeLine
//                            (
//                                dual,
//                                nearestPointOnTet
//                            );
//                            snapping2.writeLine
//                            (
//                                nearestPointOnTet,
//                                hitInfo.point()
//                            );
//
//                            pts[cit->cellIndex()] = hitInfo.point();
//                        }
//                    }
//                }
//            }

            boundaryPts[cit->cellIndex()] = classifyBoundaryPoint(cit);
        }
        else
        {
            cit->cellIndex() = Cb::ctFar;
        }
    }

    //pts.setSize(this->cellCount());

    //boundaryPts.setSize(this->cellCount());
}


void Foam::conformalVoronoiMesh::reindexDualVertices
(
    const Map<label>& dualPtIndexMap,
    labelList& boundaryPts
)
{
    for
    (
        Delaunay::Finite_cells_iterator cit = finite_cells_begin();
        cit != finite_cells_end();
        ++cit
    )
    {
        if (dualPtIndexMap.found(cit->cellIndex()))
        {
            cit->cellIndex() = dualPtIndexMap[cit->cellIndex()];
            boundaryPts[cit->cellIndex()] =
                max
                (
                    boundaryPts[cit->cellIndex()],
                    boundaryPts[dualPtIndexMap[cit->cellIndex()]]
                );
        }
    }
}


Foam::label Foam::conformalVoronoiMesh::createPatchInfo
(
    wordList& patchNames,
    PtrList<dictionary>& patchDicts
) const
{
    patchNames = geometryToConformTo_.patchNames();

    patchDicts.setSize(patchNames.size() + 1);

    const PtrList<dictionary>& patchInfo = geometryToConformTo_.patchInfo();

    forAll(patchNames, patchi)
    {
        if (patchInfo.set(patchi))
        {
            patchDicts.set(patchi, new dictionary(patchInfo[patchi]));
        }
        else
        {
            patchDicts.set(patchi, new dictionary());
            patchDicts[patchi].set
            (
                "type",
                wallPolyPatch::typeName
            );
        }
    }

    patchNames.setSize(patchNames.size() + 1);
    label defaultPatchIndex = patchNames.size() - 1;
    patchNames[defaultPatchIndex] = "foamyHexMesh_defaultPatch";
    patchDicts.set(defaultPatchIndex, new dictionary());
    patchDicts[defaultPatchIndex].set
    (
        "type",
        wallPolyPatch::typeName
    );

    label nProcPatches = 0;

    if (Pstream::parRun())
    {
        List<boolList> procUsedList
        (
            Pstream::nProcs(),
            boolList(Pstream::nProcs(), false)
        );

        boolList& procUsed = procUsedList[Pstream::myProcNo()];

        // Determine which processor patches are required
        for
        (
            Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
            vit != finite_vertices_end();
            vit++
        )
        {
            // This test is not sufficient if one of the processors does
            // not receive a referred vertex from another processor, but does
            // send one to the other processor.
            if (vit->referred())
            {
                procUsed[vit->procIndex()] = true;
            }
        }

        // Because the previous test was insufficient, combine the lists.
        Pstream::allGatherList(procUsedList);

        forAll(procUsedList, proci)
        {
            if (proci != Pstream::myProcNo())
            {
                if (procUsedList[proci][Pstream::myProcNo()])
                {
                    procUsed[proci] = true;
                }
            }
        }

        forAll(procUsed, pUI)
        {
            if (procUsed[pUI])
            {
                nProcPatches++;
            }
        }

        label nNonProcPatches = patchNames.size();
        label nTotalPatches = nNonProcPatches + nProcPatches;

        patchNames.setSize(nTotalPatches);
        patchDicts.setSize(nTotalPatches);
        for (label pI = nNonProcPatches; pI < nTotalPatches; ++pI)
        {
            patchDicts.set(pI, new dictionary());
        }

        label procAddI = 0;

        forAll(procUsed, pUI)
        {
            if (procUsed[pUI])
            {
                patchNames[nNonProcPatches + procAddI] =
                    processorPolyPatch::newName(Pstream::myProcNo(), pUI);

                patchDicts[nNonProcPatches + procAddI].set
                (
                    "type",
                    processorPolyPatch::typeName
                );

                patchDicts[nNonProcPatches + procAddI].set
                (
                    "myProcNo",
                    Pstream::myProcNo()
                );

                patchDicts[nNonProcPatches + procAddI].set("neighbProcNo", pUI);

                procAddI++;
            }
        }
    }

    return defaultPatchIndex;
}


Foam::vector Foam::conformalVoronoiMesh::calcSharedPatchNormal
(
    Cell_handle c1,
    Cell_handle c2
) const
{
    List<Foam::point> patchEdge(2, point::max);

    // Get shared Facet
    for (label cI = 0; cI < 4; ++cI)
    {
        if (c1->neighbor(cI) != c2 && !c1->vertex(cI)->constrained())
        {
            if (c1->vertex(cI)->internalBoundaryPoint())
            {
                patchEdge[0] = topoint(c1->vertex(cI)->point());
            }
            else
            {
                patchEdge[1] = topoint(c1->vertex(cI)->point());
            }
        }
    }

    Info<< "    " << patchEdge << endl;

    return vector(patchEdge[1] - patchEdge[0]);
}


bool Foam::conformalVoronoiMesh::boundaryDualFace
(
    Cell_handle c1,
    Cell_handle c2
) const
{
    label nInternal = 0;
    label nExternal = 0;

    for (label cI = 0; cI < 4; ++cI)
    {
        if (c1->neighbor(cI) != c2 && !c1->vertex(cI)->constrained())
        {
            if (c1->vertex(cI)->internalBoundaryPoint())
            {
                nInternal++;
            }
            else if (c1->vertex(cI)->externalBoundaryPoint())
            {
                nExternal++;
            }
        }
    }

    Info<< "in = " << nInternal << " out = " << nExternal << endl;

    return (nInternal == 1 && nExternal == 1);
}


void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
(
    const pointField& pts,
    faceList& faces,
    labelList& owner,
    labelList& neighbour,
    wordList& patchNames,
    PtrList<dictionary>& patchDicts,
    labelListList& patchPointPairSlaves,
    bitSet& boundaryFacesToRemove,
    bool includeEmptyPatches
) const
{
    const label defaultPatchIndex = createPatchInfo(patchNames, patchDicts);

    const label nPatches = patchNames.size();

    labelList procNeighbours(nPatches);
    forAll(procNeighbours, patchi)
    {
        procNeighbours[patchi] =
            patchDicts[patchi].getOrDefault<label>("neighbProcNo", -1);
    }

    List<DynamicList<face>> patchFaces(nPatches);
    List<DynamicList<label>> patchOwners(nPatches);
    // Per patch face the index of the slave node of the point pair
    List<DynamicList<label>> patchPPSlaves(nPatches);
    List<DynamicList<bool>> indirectPatchFace(nPatches);


    faces.setSize(number_of_finite_edges());
    owner.setSize(number_of_finite_edges());
    neighbour.setSize(number_of_finite_edges());
    boundaryFacesToRemove.setSize(number_of_finite_edges(), false);

    labelPairPairDynListList procPatchSortingIndex(nPatches);

    label dualFacei = 0;

    if (foamyHexMeshControls().guardFeaturePoints())
    {
        OBJstream startCellStr("startingCell.obj");
        OBJstream featurePointFacesStr("ftPtFaces.obj");
        OBJstream featurePointDualsStr("ftPtDuals.obj");
        OFstream cellStr("vertexCells.obj");

        label vcount = 1;

        for
        (
            Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
            vit != finite_vertices_end();
            ++vit
        )
        {
            if (vit->constrained())
            {
                // Find a starting cell
                std::list<Cell_handle> vertexCells;
                finite_incident_cells(vit, std::back_inserter(vertexCells));

                Cell_handle startCell;

                for
                (
                    std::list<Cell_handle>::iterator vcit = vertexCells.begin();
                    vcit != vertexCells.end();
                    ++vcit
                )
                {
                    if ((*vcit)->featurePointExternalCell())
                    {
                        startCell = *vcit;
                    }

                    if ((*vcit)->real())
                    {
                        featurePointDualsStr.writeLine
                        (
                            topoint(vit->point()),
                            (*vcit)->dual()
                        );
                    }
                }

                // Error if startCell is null
                if (startCell == nullptr)
                {
                    Pout<< "Start cell is null!" << endl;
                }

                // Need to pick a direction to walk in
                Cell_handle vc1 = startCell;
                Cell_handle vc2;

                Info<< "c1 index = " << vc1->cellIndex() << " "
                    << vc1->dual() << endl;

                for (label cI = 0; cI < 4; ++cI)
                {
                    Info<< "c1 = " << cI << " "
                        << vc1->neighbor(cI)->cellIndex() << " v = "
                        << vc1->neighbor(cI)->dual() << endl;

                    Info<< vc1->vertex(cI)->info();
                }

                Cell_handle nextCell;

                for (label cI = 0; cI < 4; ++cI)
                {
                    if (vc1->vertex(cI)->externalBoundaryPoint())
                    {
                        vc2 = vc1->neighbor(cI);

                        Info<< "    c2 is neighbor "
                            << vc2->cellIndex()
                            << " of c1" << endl;

                        for (label cI = 0; cI < 4; ++cI)
                        {
                            Info<< "    c2 = " << cI << " "
                                << vc2->neighbor(cI)->cellIndex() << " v = "
                                << vc2->vertex(cI)->index() << endl;
                        }

                        face f(3);
                        f[0] = vit->index();
                        f[1] = vc1->cellIndex();
                        f[2] = vc2->cellIndex();

                        Info<< "f " << f << endl;
                        forAll(f, pI)
                        {
                            Info<< "    " << pts[f[pI]] << endl;
                        }

                        vector correctNormal = calcSharedPatchNormal(vc1, vc2);
                        correctNormal.normalise();

                        Info<< "    cN " << correctNormal << endl;

                        vector fN = f.areaNormal(pts);

                        if (mag(fN) < SMALL)
                        {
                            nextCell = vc2;
                            continue;
                        }

                        fN.normalise();
                        Info<< "    fN " << fN << endl;

                        if ((fN & correctNormal) > 0)
                        {
                            nextCell = vc2;
                            break;
                        }
                    }
                }

                vc2 = nextCell;

                label own = vit->index();
                face f(3);
                f[0] = own;

                Info<< "Start walk from " << vc1->cellIndex()
                    << " to " << vc2->cellIndex() << endl;

                // Walk while not at start cell

                label iter = 0;
                do
                {
                    Info<< "     Walk from " << vc1->cellIndex()
                        << " " << vc1->dual()
                        << " to " << vc2->cellIndex()
                        << " " << vc2->dual()
                        << endl;

                    startCellStr.writeLine(vc1->dual(), vc2->dual());

                    // Get patch by getting face between cells and the two
                    // points on the face that are not the feature vertex
                    label patchIndex =
                        geometryToConformTo_.findPatch
                        (
                            topoint(vit->point())
                        );

                    f[1] = vc1->cellIndex();
                    f[2] = vc2->cellIndex();

                    patchFaces[patchIndex].append(f);
                    patchOwners[patchIndex].append(own);
                    patchPPSlaves[patchIndex].append(own);

                    // Find next cell
                    Cell_handle nextCell;

                    Info<< "    c1 vertices " << vc2->dual() << endl;
                    for (label cI = 0; cI < 4; ++cI)
                    {
                        Info<< "        " << vc2->vertex(cI)->info();
                    }
                    Info<< "    c1 neighbour vertices " << endl;
                    for (label cI = 0; cI < 4; ++cI)
                    {
                        if
                        (
                            !vc2->vertex(cI)->constrained()
                         && vc2->neighbor(cI) != vc1
                         && !is_infinite(vc2->neighbor(cI))
                         &&
                            (
                                vc2->neighbor(cI)->featurePointExternalCell()
                             || vc2->neighbor(cI)->featurePointInternalCell()
                            )
                         && vc2->neighbor(cI)->hasConstrainedPoint()
                        )
                        {
                            DelaunayMeshTools::drawDelaunayCell
                            (
                                cellStr,
                                vc2->neighbor(cI),
                                vcount++
                            );

                            Info<< "        neighbour " << cI << " "
                                << vc2->neighbor(cI)->dual() << endl;
                            for (label I = 0; I < 4; ++I)
                            {
                                Info<< "            "
                                    << vc2->neighbor(cI)->vertex(I)->info();
                            }
                        }
                    }

                    for (label cI = 0; cI < 4; ++cI)
                    {
                        if
                        (
                            !vc2->vertex(cI)->constrained()
                         && vc2->neighbor(cI) != vc1
                         && !is_infinite(vc2->neighbor(cI))
                         &&
                            (
                                vc2->neighbor(cI)->featurePointExternalCell()
                             || vc2->neighbor(cI)->featurePointInternalCell()
                            )
                         && vc2->neighbor(cI)->hasConstrainedPoint()
                        )
                        {
                            // check if shared edge is internal/internal
                            if (boundaryDualFace(vc2, vc2->neighbor(cI)))
                            {
                                nextCell = vc2->neighbor(cI);
                                break;
                            }
                        }
                    }

                    vc1 = vc2;
                    vc2 = nextCell;

                    iter++;
                } while (vc1 != startCell && iter < 100);
            }
        }
    }

    for
    (
        Delaunay::Finite_edges_iterator eit = finite_edges_begin();
        eit != finite_edges_end();
        ++eit
    )
    {
        Cell_handle c = eit->first;
        Vertex_handle vA = c->vertex(eit->second);
        Vertex_handle vB = c->vertex(eit->third);

        if (vA->constrained() && vB->constrained())
        {
            continue;
        }

        if
        (
            (vA->constrained() && vB->internalOrBoundaryPoint())
         || (vB->constrained() && vA->internalOrBoundaryPoint())
        )
        {
            face newDualFace = buildDualFace(eit);

            label own = -1;
            label nei = -1;

            if (ownerAndNeighbour(vA, vB, own, nei))
            {
                reverse(newDualFace);
            }

            // internal face
            faces[dualFacei] = newDualFace;
            owner[dualFacei] = own;
            neighbour[dualFacei] = nei;

            dualFacei++;
        }
        else if
        (
            (vA->internalOrBoundaryPoint() && !vA->referred())
         || (vB->internalOrBoundaryPoint() && !vB->referred())
        )
        {
            if
            (
                (vA->internalPoint() && vB->externalBoundaryPoint())
             || (vB->internalPoint() && vA->externalBoundaryPoint())
            )
            {
                Cell_circulator ccStart = incident_cells(*eit);
                Cell_circulator cc1 = ccStart;
                Cell_circulator cc2 = cc1;

                cc2++;

                bool skipEdge = false;

                do
                {
                    if
                    (
                        cc1->hasFarPoint() || cc2->hasFarPoint()
                     || is_infinite(cc1) || is_infinite(cc2)
                    )
                    {
                        Pout<< "Ignoring edge between internal and external: "
                            << vA->info()
                            << vB->info();

                        skipEdge = true;
                        break;
                    }

                    cc1++;
                    cc2++;

                } while (cc1 != ccStart);


                // Do not create faces if the internal point is outside!
                // This occurs because the internal point is not determined to
                // be outside in the inside/outside test. This is most likely
                // due to the triangle.nearestPointClassify test not returning
                // edge/point as the nearest type.

                if (skipEdge)
                {
                    continue;
                }
            }

            face newDualFace = buildDualFace(eit);

            if (newDualFace.size() >= 3)
            {
                label own = -1;
                label nei = -1;

                if (ownerAndNeighbour(vA, vB, own, nei))
                {
                    reverse(newDualFace);
                }

                label patchIndex = -1;

                pointFromPoint ptA = topoint(vA->point());
                pointFromPoint ptB = topoint(vB->point());

                if (nei == -1)
                {
                    // boundary face

                    if (isProcBoundaryEdge(eit))
                    {
                        // One (and only one) of the points is an internal
                        // point from another processor

                        label procIndex = max(vA->procIndex(), vB->procIndex());

                        patchIndex = max
                        (
                            procNeighbours.find(vA->procIndex()),
                            procNeighbours.find(vB->procIndex())
                        );

                        // The lower processor index is the owner of the
                        // two for the purpose of sorting the patch faces.

                        if (Pstream::myProcNo() < procIndex)
                        {
                            // Use this processor's vertex index as the master
                            // for sorting

                            DynamicList<labelPairPair>& sortingIndex =
                                procPatchSortingIndex[patchIndex];

                            if (vB->internalOrBoundaryPoint() && vB->referred())
                            {
                                sortingIndex.append
                                (
                                    labelPairPair
                                    (
                                        labelPair(vA->index(), vA->procIndex()),
                                        labelPair(vB->index(), vB->procIndex())
                                    )
                                );
                            }
                            else
                            {
                                sortingIndex.append
                                (
                                    labelPairPair
                                    (
                                        labelPair(vB->index(), vB->procIndex()),
                                        labelPair(vA->index(), vA->procIndex())
                                    )
                                );
                            }
                        }
                        else
                        {
                            // Use the other processor's vertex index as the
                            // master for sorting

                            DynamicList<labelPairPair>& sortingIndex =
                                procPatchSortingIndex[patchIndex];

                            if (vA->internalOrBoundaryPoint() && vA->referred())
                            {
                                sortingIndex.append
                                (
                                    labelPairPair
                                    (
                                        labelPair(vA->index(), vA->procIndex()),
                                        labelPair(vB->index(), vB->procIndex())
                                    )
                                );
                            }
                            else
                            {
                                sortingIndex.append
                                (
                                    labelPairPair
                                    (
                                        labelPair(vB->index(), vB->procIndex()),
                                        labelPair(vA->index(), vA->procIndex())
                                    )
                                );
                            }
                        }

//                        Pout<< ptA << " " << ptB
//                            << " proc indices "
//                            << vA->procIndex() << " " << vB->procIndex()
//                            << " indices " << vA->index()
//                            << " " << vB->index()
//                            << " my proc " << Pstream::myProcNo()
//                            << " addedIndex "
//                            << procPatchSortingIndex[patchIndex].last()
//                            << endl;
                    }
                    else
                    {
                        patchIndex = geometryToConformTo_.findPatch(ptA, ptB);
                    }

                    if (patchIndex == -1)
                    {
                        // Did not find a surface patch between
                        // between Dv pair, finding nearest patch

//                         Pout<< "Did not find a surface patch between "
//                             << "for face, finding nearest patch to"
//                             << 0.5*(ptA + ptB) << endl;

                        patchIndex = geometryToConformTo_.findPatch
                        (
                            0.5*(ptA + ptB)
                        );
                    }

                    patchFaces[patchIndex].append(newDualFace);
                    patchOwners[patchIndex].append(own);

                    // If the two vertices are a pair, then the patch face is
                    // a desired one.
                    if
                    (
                        vA->boundaryPoint() && vB->boundaryPoint()
                     && !ptPairs_.isPointPair(vA, vB)
                     && !ftPtConformer_.featurePointPairs().isPointPair(vA, vB)
                    )
                    {
                        indirectPatchFace[patchIndex].append(true);
                    }
                    else
                    {
                        indirectPatchFace[patchIndex].append(false);
                    }

                    // Store the non-internal or boundary point
                    if (vA->internalOrBoundaryPoint())
                    {
                        patchPPSlaves[patchIndex].append(vB->index());
                    }
                    else
                    {
                        patchPPSlaves[patchIndex].append(vA->index());
                    }
                }
                else
                {
                    if
                    (
                        !vA->boundaryPoint()
                     || !vB->boundaryPoint()
                     || ptPairs_.isPointPair(vA, vB)
                     || ftPtConformer_.featurePointPairs().isPointPair(vA, vB)
                    )
                    {
                        patchIndex = geometryToConformTo_.findPatch(ptA, ptB);
                    }

                    if
                    (
                        patchIndex != -1
                     && geometryToConformTo_.patchInfo().set(patchIndex)
                    )
                    {
                        // baffle faces

                        patchFaces[patchIndex].append(newDualFace);
                        patchOwners[patchIndex].append(own);
                        indirectPatchFace[patchIndex].append(false);

                        reverse(newDualFace);

                        patchFaces[patchIndex].append(newDualFace);
                        patchOwners[patchIndex].append(nei);
                        indirectPatchFace[patchIndex].append(false);

                        if
                        (
                            labelPair(vB->index(), vB->procIndex())
                          < labelPair(vA->index(), vA->procIndex())
                        )
                        {
                            patchPPSlaves[patchIndex].append(vB->index());
                            patchPPSlaves[patchIndex].append(vB->index());
                        }
                        else
                        {
                            patchPPSlaves[patchIndex].append(vA->index());
                            patchPPSlaves[patchIndex].append(vA->index());
                        }

                    }
                    else
                    {
                        // internal face
                        faces[dualFacei] = newDualFace;
                        owner[dualFacei] = own;
                        neighbour[dualFacei] = nei;

                        dualFacei++;
                    }
                }
            }
        }
    }

    if (!patchFaces[defaultPatchIndex].empty())
    {
        Pout<< nl << patchFaces[defaultPatchIndex].size()
            << " faces were not able to have their patch determined from "
            << "the surface. "
            << nl <<  "Adding to patch " << patchNames[defaultPatchIndex]
            << endl;
    }

    label nInternalFaces = dualFacei;

    faces.setSize(nInternalFaces);
    owner.setSize(nInternalFaces);
    neighbour.setSize(nInternalFaces);

    timeCheck("polyMesh quality checked");

    sortFaces(faces, owner, neighbour);

    sortProcPatches
    (
        patchFaces,
        patchOwners,
        patchPPSlaves,
        procPatchSortingIndex
    );

    timeCheck("faces, owner, neighbour sorted");

    addPatches
    (
        nInternalFaces,
        faces,
        owner,
        patchDicts,
        boundaryFacesToRemove,
        patchFaces,
        patchOwners,
        indirectPatchFace
    );

    // Return     patchPointPairSlaves.setSize(nPatches);
    patchPointPairSlaves.setSize(nPatches);
    forAll(patchPPSlaves, patchi)
    {
        patchPointPairSlaves[patchi].transfer(patchPPSlaves[patchi]);
    }

    if (foamyHexMeshControls().objOutput())
    {
        Info<< "Writing processor interfaces" << endl;

        forAll(patchDicts, nbI)
        {
            if (patchFaces[nbI].size() > 0)
            {
                const label neighbour =
                    patchDicts[nbI].getOrDefault<label>("neighbProcNo", -1);

                faceList procPatchFaces = patchFaces[nbI];

                // Reverse faces as it makes it easier to analyse the output
                // using a diff
                if (neighbour < Pstream::myProcNo())
                {
                    forAll(procPatchFaces, fI)
                    {
                        procPatchFaces[fI].flip();
                    }
                }

                if (neighbour != -1)
                {
                    word fName =
                        "processor_"
                      + name(Pstream::myProcNo())
                      + "_to_"
                      + name(neighbour)
                      + "_interface.obj";

                    DelaunayMeshTools::writeProcessorInterface
                    (
                        time().path()/fName,
                        *this,
                        procPatchFaces
                    );
                }
            }
        }
    }
}


void Foam::conformalVoronoiMesh::sortFaces
(
    faceList& faces,
    labelList& owner,
    labelList& neighbour
) const
{
    // Upper triangular order:
    // + owner is sorted in ascending cell order
    // + within each block of equal value for owner, neighbour is sorted in
    //   ascending cell order.
    // + faces sorted to correspond
    // e.g.
    // owner | neighbour
    // 0     | 2
    // 0     | 23
    // 0     | 71
    // 1     | 23
    // 1     | 24
    // 1     | 91

    List<labelPair> ownerNeighbourPair(owner.size());

    forAll(ownerNeighbourPair, oNI)
    {
        ownerNeighbourPair[oNI] = labelPair(owner[oNI], neighbour[oNI]);
    }

    Info<< nl
        << "Sorting faces, owner and neighbour into upper triangular order"
        << endl;

    labelList oldToNew(sortedOrder(ownerNeighbourPair));
    oldToNew = invert(oldToNew.size(), oldToNew);

    inplaceReorder(oldToNew, faces);
    inplaceReorder(oldToNew, owner);
    inplaceReorder(oldToNew, neighbour);
}


void Foam::conformalVoronoiMesh::sortProcPatches
(
    List<DynamicList<face>>& patchFaces,
    List<DynamicList<label>>& patchOwners,
    List<DynamicList<label>>& patchPointPairSlaves,
    labelPairPairDynListList& patchSortingIndices
) const
{
    if (!Pstream::parRun())
    {
        return;
    }

    forAll(patchSortingIndices, patchi)
    {
        faceList& faces = patchFaces[patchi];
        labelList& owner = patchOwners[patchi];
        DynamicList<label>& slaves = patchPointPairSlaves[patchi];
        DynamicList<labelPairPair>& sortingIndices
            = patchSortingIndices[patchi];

        if (!sortingIndices.empty())
        {
            if
            (
                faces.size() != sortingIndices.size()
             || owner.size() != sortingIndices.size()
             || slaves.size() != sortingIndices.size()
            )
            {
                FatalErrorInFunction
                    << "patch size and size of sorting indices is inconsistent "
                    << " for patch " << patchi << nl
                    << " faces.size() " << faces.size() << nl
                    << " owner.size() " << owner.size() << nl
                    << " slaves.size() " << slaves.size() << nl
                    << " sortingIndices.size() "
                    << sortingIndices.size()
                    << exit(FatalError) << endl;
            }

            labelList oldToNew(sortedOrder(sortingIndices));
            oldToNew = invert(oldToNew.size(), oldToNew);

            inplaceReorder(oldToNew, sortingIndices);
            inplaceReorder(oldToNew, faces);
            inplaceReorder(oldToNew, owner);
            inplaceReorder(oldToNew, slaves);
        }
    }
}


void Foam::conformalVoronoiMesh::addPatches
(
    const label nInternalFaces,
    faceList& faces,
    labelList& owner,
    PtrList<dictionary>& patchDicts,
    bitSet& boundaryFacesToRemove,
    const List<DynamicList<face>>& patchFaces,
    const List<DynamicList<label>>& patchOwners,
    const List<DynamicList<bool>>& indirectPatchFace
) const
{
    label nBoundaryFaces = 0;

    forAll(patchFaces, p)
    {
        patchDicts[p].set("nFaces", patchFaces[p].size());
        patchDicts[p].set("startFace", nInternalFaces + nBoundaryFaces);

        nBoundaryFaces += patchFaces[p].size();
    }

    faces.setSize(nInternalFaces + nBoundaryFaces);
    owner.setSize(nInternalFaces + nBoundaryFaces);
    boundaryFacesToRemove.setSize(nInternalFaces + nBoundaryFaces);

    label facei = nInternalFaces;

    forAll(patchFaces, p)
    {
        forAll(patchFaces[p], f)
        {
            faces[facei] = patchFaces[p][f];
            owner[facei] = patchOwners[p][f];
            boundaryFacesToRemove[facei] = indirectPatchFace[p][f];

            facei++;
        }
    }
}


void Foam::conformalVoronoiMesh::removeUnusedPoints
(
    faceList& faces,
    pointField& pts,
    labelList& boundaryPts
) const
{
    Info<< nl << "Removing unused points" << endl;

    bitSet ptUsed(pts.size(), false);

    // Scan all faces to find all of the points that are used

    forAll(faces, fI)
    {
        const face& f = faces[fI];

        ptUsed.set(f);
    }

    label pointi = 0;

    labelList oldToNew(pts.size(), label(-1));

    // Move all of the used points to the start of the pointField and
    // truncate it

    forAll(ptUsed, ptUI)
    {
        if (ptUsed.test(ptUI))
        {
            oldToNew[ptUI] = pointi++;
        }
    }

    inplaceReorder(oldToNew, pts);
    inplaceReorder(oldToNew, boundaryPts);

    Info<< "    Removing "
        << returnReduce(pts.size() - pointi, sumOp<label>())
        << " unused points"
        << endl;

    pts.setSize(pointi);
    boundaryPts.setSize(pointi);

    // Renumber the faces to use the new point numbers

    forAll(faces, fI)
    {
        inplaceRenumber(oldToNew, faces[fI]);
    }
}


Foam::labelList Foam::conformalVoronoiMesh::removeUnusedCells
(
    labelList& owner,
    labelList& neighbour
) const
{
    Info<< nl << "Removing unused cells" << endl;

    bitSet cellUsed(vertexCount(), false);

    // Scan all faces to find all of the cells that are used

    cellUsed.set(owner);
    cellUsed.set(neighbour);

    label celli = 0;

    labelList oldToNew(cellUsed.size(), label(-1));

    // Move all of the used cellCentres to the start of the pointField and
    // truncate it

    forAll(cellUsed, cellUI)
    {
        if (cellUsed.test(cellUI))
        {
            oldToNew[cellUI] = celli++;
        }
    }

    labelList newToOld(invert(celli, oldToNew));

    // Find all of the unused cells, create a list of them, then
    // subtract one from each owner and neighbour entry for each of
    // the unused cell indices that it is above.

    DynamicList<label> unusedCells;

    forAll(cellUsed, cUI)
    {
        if (!cellUsed.test(cUI))
        {
            unusedCells.append(cUI);
        }
    }

    if (unusedCells.size() > 0)
    {
        Info<< "    Removing "
            << returnReduce(unusedCells.size(), sumOp<label>())
            <<  " unused cell labels" << endl;

        forAll(owner, oI)
        {
            label& o = owner[oI];

            o -= findLower(unusedCells, o) + 1;
        }

        forAll(neighbour, nI)
        {
            label& n = neighbour[nI];

            n -= findLower(unusedCells, n) + 1;
        }
    }

    return newToOld;
}


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