autoLayerDriver.C

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

Description
    All to do with adding cell layers

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

#include "autoLayerDriver.H"
#include "fvMesh.H"
#include "Time.H"
#include "meshRefinement.H"
#include "removePoints.H"
#include "pointFields.H"
#include "motionSmoother.H"
#include "unitConversion.H"
#include "pointSet.H"
#include "faceSet.H"
#include "cellSet.H"
#include "polyTopoChange.H"
#include "mapPolyMesh.H"
#include "addPatchCellLayer.H"
#include "mapDistributePolyMesh.H"
#include "OBJstream.H"
#include "layerParameters.H"
#include "combineFaces.H"
#include "IOmanip.H"
#include "globalIndex.H"
#include "DynamicField.H"
#include "PatchTools.H"
#include "slipPointPatchFields.H"
#include "fixedValuePointPatchFields.H"
#include "zeroFixedValuePointPatchFields.H"
#include "calculatedPointPatchFields.H"
#include "cyclicSlipPointPatchFields.H"
#include "fixedValueFvPatchFields.H"
#include "localPointRegion.H"
#include "externalDisplacementMeshMover.H"
#include "scalarIOField.H"

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

namespace Foam
{

defineTypeNameAndDebug(autoLayerDriver, 0);

} // End namespace Foam


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

// For debugging: Dump displacement to .obj files
void Foam::autoLayerDriver::dumpDisplacement
(
    const fileName& prefix,
    const indirectPrimitivePatch& pp,
    const vectorField& patchDisp,
    const List<extrudeMode>& extrudeStatus
)
{
    OBJstream dispStr(prefix + "_disp.obj");
    Info<< "Writing all displacements to " << dispStr.name() << endl;

    forAll(patchDisp, patchPointI)
    {
        const point& pt = pp.localPoints()[patchPointI];
        dispStr.write(linePointRef(pt, pt + patchDisp[patchPointI]));
    }


    OBJstream illStr(prefix + "_illegal.obj");
    Info<< "Writing invalid displacements to " << illStr.name() << endl;

    forAll(patchDisp, patchPointI)
    {
        if (extrudeStatus[patchPointI] != EXTRUDE)
        {
            const point& pt = pp.localPoints()[patchPointI];
            illStr.write(linePointRef(pt, pt + patchDisp[patchPointI]));
        }
    }
}


Foam::tmp<Foam::scalarField> Foam::autoLayerDriver::avgPointData
(
    const indirectPrimitivePatch& pp,
    const scalarField& pointFld
)
{
    tmp<scalarField> tfaceFld(new scalarField(pp.size(), 0.0));
    scalarField& faceFld = tfaceFld();

    forAll(pp.localFaces(), faceI)
    {
        const face& f = pp.localFaces()[faceI];
        if (f.size())
        {
            forAll(f, fp)
            {
                faceFld[faceI] += pointFld[f[fp]];
            }
            faceFld[faceI] /= f.size();
        }
    }
    return tfaceFld;
}


// Check that primitivePatch is not multiply connected. Collect non-manifold
// points in pointSet.
void Foam::autoLayerDriver::checkManifold
(
    const indirectPrimitivePatch& fp,
    pointSet& nonManifoldPoints
)
{
    // Check for non-manifold points (surface pinched at point)
    fp.checkPointManifold(false, &nonManifoldPoints);

    // Check for edge-faces (surface pinched at edge)
    const labelListList& edgeFaces = fp.edgeFaces();

    forAll(edgeFaces, edgeI)
    {
        const labelList& eFaces = edgeFaces[edgeI];

        if (eFaces.size() > 2)
        {
            const edge& e = fp.edges()[edgeI];

            nonManifoldPoints.insert(fp.meshPoints()[e[0]]);
            nonManifoldPoints.insert(fp.meshPoints()[e[1]]);
        }
    }
}


void Foam::autoLayerDriver::checkMeshManifold() const
{
    const fvMesh& mesh = meshRefiner_.mesh();

    Info<< nl << "Checking mesh manifoldness ..." << endl;

    // Get all outside faces
    labelList outsideFaces(mesh.nFaces() - mesh.nInternalFaces());

    for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
    {
         outsideFaces[faceI - mesh.nInternalFaces()] = faceI;
    }

    pointSet nonManifoldPoints
    (
        mesh,
        "nonManifoldPoints",
        mesh.nPoints() / 100
    );

    // Build primitivePatch out of faces and check it for problems.
    checkManifold
    (
        indirectPrimitivePatch
        (
            IndirectList<face>(mesh.faces(), outsideFaces),
            mesh.points()
        ),
        nonManifoldPoints
    );

    label nNonManif = returnReduce(nonManifoldPoints.size(), sumOp<label>());

    if (nNonManif > 0)
    {
        Info<< "Outside of mesh is multiply connected across edges or"
            << " points." << nl
            << "This is not a fatal error but might cause some unexpected"
            << " behaviour." << nl
            //<< "Writing " << nNonManif
            //<< " points where this happens to pointSet "
            //<< nonManifoldPoints.name()
            << endl;

        //nonManifoldPoints.instance() = meshRefiner_.timeName();
        //nonManifoldPoints.write();
    }
    Info<< endl;
}



// Unset extrusion on point. Returns true if anything unset.
bool Foam::autoLayerDriver::unmarkExtrusion
(
    const label patchPointI,
    pointField& patchDisp,
    labelList& patchNLayers,
    List<extrudeMode>& extrudeStatus
)
{
    if (extrudeStatus[patchPointI] == EXTRUDE)
    {
        extrudeStatus[patchPointI] = NOEXTRUDE;
        patchNLayers[patchPointI] = 0;
        patchDisp[patchPointI] = vector::zero;
        return true;
    }
    else if (extrudeStatus[patchPointI] == EXTRUDEREMOVE)
    {
        extrudeStatus[patchPointI] = NOEXTRUDE;
        patchNLayers[patchPointI] = 0;
        patchDisp[patchPointI] = vector::zero;
        return true;
    }
    else
    {
        return false;
    }
}


// Unset extrusion on face. Returns true if anything unset.
bool Foam::autoLayerDriver::unmarkExtrusion
(
    const face& localFace,
    pointField& patchDisp,
    labelList& patchNLayers,
    List<extrudeMode>& extrudeStatus
)
{
    bool unextruded = false;

    forAll(localFace, fp)
    {
        if
        (
            unmarkExtrusion
            (
                localFace[fp],
                patchDisp,
                patchNLayers,
                extrudeStatus
            )
        )
        {
            unextruded = true;
        }
    }
    return unextruded;
}


// No extrusion at non-manifold points.
void Foam::autoLayerDriver::handleNonManifolds
(
    const indirectPrimitivePatch& pp,
    const labelList& meshEdges,
    const labelListList& edgeGlobalFaces,
    pointField& patchDisp,
    labelList& patchNLayers,
    List<extrudeMode>& extrudeStatus
) const
{
    const fvMesh& mesh = meshRefiner_.mesh();

    Info<< nl << "Handling non-manifold points ..." << endl;

    // Detect non-manifold points
    Info<< nl << "Checking patch manifoldness ..." << endl;

    pointSet nonManifoldPoints(mesh, "nonManifoldPoints", pp.nPoints());

    // 1. Local check
    checkManifold(pp, nonManifoldPoints);

    // 2. Remote check for boundary edges on coupled boundaries
    forAll(edgeGlobalFaces, edgeI)
    {
        if
        (
            pp.edgeFaces()[edgeI].size() == 1
         && edgeGlobalFaces[edgeI].size() > 2
        )
        {
            // So boundary edges that are connected to more than 2 processors
            // i.e. a non-manifold edge which is exactly on a processor
            // boundary.
            const edge& e = pp.edges()[edgeI];
            nonManifoldPoints.insert(pp.meshPoints()[e[0]]);
            nonManifoldPoints.insert(pp.meshPoints()[e[1]]);
        }
    }

    // 3. Remote check for end of layer across coupled boundaries
    {
        PackedBoolList isCoupledEdge(mesh.nEdges());

        const labelList& cpEdges = mesh.globalData().coupledPatchMeshEdges();
        forAll(cpEdges, i)
        {
            isCoupledEdge[cpEdges[i]] = true;
        }
        syncTools::syncEdgeList
        (
            mesh,
            isCoupledEdge,
            orEqOp<unsigned int>(),
            0
        );

        forAll(edgeGlobalFaces, edgeI)
        {
            label meshEdgeI = meshEdges[edgeI];

            if
            (
                pp.edgeFaces()[edgeI].size() == 1
             && edgeGlobalFaces[edgeI].size() == 1
             && isCoupledEdge[meshEdgeI]
            )
            {
                // Edge of patch but no continuation across processor.
                const edge& e = pp.edges()[edgeI];
                //Pout<< "** Stopping extrusion on edge "
                //    << pp.localPoints()[e[0]]
                //    << pp.localPoints()[e[1]] << endl;
                nonManifoldPoints.insert(pp.meshPoints()[e[0]]);
                nonManifoldPoints.insert(pp.meshPoints()[e[1]]);
            }
        }
    }



    label nNonManif = returnReduce(nonManifoldPoints.size(), sumOp<label>());

    Info<< "Outside of local patch is multiply connected across edges or"
        << " points at " << nNonManif << " points." << endl;

    if (nNonManif > 0)
    {
        const labelList& meshPoints = pp.meshPoints();

        forAll(meshPoints, patchPointI)
        {
            if (nonManifoldPoints.found(meshPoints[patchPointI]))
            {
                unmarkExtrusion
                (
                    patchPointI,
                    patchDisp,
                    patchNLayers,
                    extrudeStatus
                );
            }
        }
    }

    Info<< "Set displacement to zero for all " << nNonManif
        << " non-manifold points" << endl;
}


// Parallel feature edge detection. Assumes non-manifold edges already handled.
void Foam::autoLayerDriver::handleFeatureAngle
(
    const indirectPrimitivePatch& pp,
    const labelList& meshEdges,
    const scalar minCos,
    pointField& patchDisp,
    labelList& patchNLayers,
    List<extrudeMode>& extrudeStatus
) const
{
    const fvMesh& mesh = meshRefiner_.mesh();

    Info<< nl << "Handling feature edges ..." << endl;

    if (minCos < 1-SMALL)
    {
        // Normal component of normals of connected faces.
        vectorField edgeNormal(mesh.nEdges(), point::max);

        const labelListList& edgeFaces = pp.edgeFaces();

        forAll(edgeFaces, edgeI)
        {
            const labelList& eFaces = pp.edgeFaces()[edgeI];

            label meshEdgeI = meshEdges[edgeI];

            forAll(eFaces, i)
            {
                nomalsCombine()
                (
                    edgeNormal[meshEdgeI],
                    pp.faceNormals()[eFaces[i]]
                );
            }
        }

        syncTools::syncEdgeList
        (
            mesh,
            edgeNormal,
            nomalsCombine(),
            point::max          // null value
        );

        autoPtr<OBJstream> str;
        if (debug&meshRefinement::MESH)
        {
            str.reset
            (
                new OBJstream
                (
                    mesh.time().path()
                  / "featureEdges_"
                  + meshRefiner_.timeName()
                  + ".obj"
                )
            );
            Info<< "Writing feature edges to " << str().name() << endl;
        }

        label nFeats = 0;

        // Now on coupled edges the edgeNormal will have been truncated and
        // only be still be the old value where two faces have the same normal
        forAll(edgeFaces, edgeI)
        {
            const labelList& eFaces = pp.edgeFaces()[edgeI];

            label meshEdgeI = meshEdges[edgeI];

            const vector& n = edgeNormal[meshEdgeI];

            if (n != point::max)
            {
                scalar cos = n & pp.faceNormals()[eFaces[0]];

                if (cos < minCos)
                {
                    const edge& e = pp.edges()[edgeI];

                    unmarkExtrusion
                    (
                        e[0],
                        patchDisp,
                        patchNLayers,
                        extrudeStatus
                    );
                    unmarkExtrusion
                    (
                        e[1],
                        patchDisp,
                        patchNLayers,
                        extrudeStatus
                    );

                    nFeats++;

                    if (str.valid())
                    {
                        const point& p0 = pp.localPoints()[e[0]];
                        const point& p1 = pp.localPoints()[e[1]];
                        str().write(linePointRef(p0, p1));
                    }
                }
            }
        }

        Info<< "Set displacement to zero for points on "
            << returnReduce(nFeats, sumOp<label>())
            << " feature edges" << endl;
    }
}


// No extrusion on cells with warped faces. Calculates the thickness of the
// layer and compares it to the space the warped face takes up. Disables
// extrusion if layer thickness is more than faceRatio of the thickness of
// the face.
void Foam::autoLayerDriver::handleWarpedFaces
(
    const indirectPrimitivePatch& pp,
    const scalar faceRatio,
    const scalar edge0Len,
    const labelList& cellLevel,
    pointField& patchDisp,
    labelList& patchNLayers,
    List<extrudeMode>& extrudeStatus
) const
{
    const fvMesh& mesh = meshRefiner_.mesh();

    Info<< nl << "Handling cells with warped patch faces ..." << nl;

    const pointField& points = mesh.points();

    label nWarpedFaces = 0;

    forAll(pp, i)
    {
        const face& f = pp[i];

        if (f.size() > 3)
        {
            label faceI = pp.addressing()[i];

            label ownLevel = cellLevel[mesh.faceOwner()[faceI]];
            scalar edgeLen = edge0Len/(1<<ownLevel);

            // Normal distance to face centre plane
            const point& fc = mesh.faceCentres()[faceI];
            const vector& fn = pp.faceNormals()[i];

            scalarField vProj(f.size());

            forAll(f, fp)
            {
                vector n = points[f[fp]] - fc;
                vProj[fp] = (n & fn);
            }

            // Get normal 'span' of face
            scalar minVal = min(vProj);
            scalar maxVal = max(vProj);

            if ((maxVal - minVal) > faceRatio * edgeLen)
            {
                if
                (
                    unmarkExtrusion
                    (
                        pp.localFaces()[i],
                        patchDisp,
                        patchNLayers,
                        extrudeStatus
                    )
                )
                {
                    nWarpedFaces++;
                }
            }
        }
    }

    Info<< "Set displacement to zero on "
        << returnReduce(nWarpedFaces, sumOp<label>())
        << " warped faces since layer would be > " << faceRatio
        << " of the size of the bounding box." << endl;
}


//// No extrusion on cells with multiple patch faces. There ususally is a reason
//// why combinePatchFaces hasn't succeeded.
//void Foam::autoLayerDriver::handleMultiplePatchFaces
//(
//    const indirectPrimitivePatch& pp,
//    pointField& patchDisp,
//    labelList& patchNLayers,
//    List<extrudeMode>& extrudeStatus
//) const
//{
//    const fvMesh& mesh = meshRefiner_.mesh();
//
//    Info<< nl << "Handling cells with multiple patch faces ..." << nl;
//
//    const labelListList& pointFaces = pp.pointFaces();
//
//    // Cells that should not get an extrusion layer
//    cellSet multiPatchCells(mesh, "multiPatchCells", pp.size());
//
//    // Detect points that use multiple faces on same cell.
//    forAll(pointFaces, patchPointI)
//    {
//        const labelList& pFaces = pointFaces[patchPointI];
//
//        labelHashSet pointCells(pFaces.size());
//
//        forAll(pFaces, i)
//        {
//            label cellI = mesh.faceOwner()[pp.addressing()[pFaces[i]]];
//
//            if (!pointCells.insert(cellI))
//            {
//                // Second or more occurrence of cell so cell has two or more
//                // pp faces connected to this point.
//                multiPatchCells.insert(cellI);
//            }
//        }
//    }
//
//    label nMultiPatchCells = returnReduce
//    (
//        multiPatchCells.size(),
//        sumOp<label>()
//    );
//
//    Info<< "Detected " << nMultiPatchCells
//        << " cells with multiple (connected) patch faces." << endl;
//
//    label nChanged = 0;
//
//    if (nMultiPatchCells > 0)
//    {
//        multiPatchCells.instance() = meshRefiner_.timeName();
//        Info<< "Writing " << nMultiPatchCells
//            << " cells with multiple (connected) patch faces to cellSet "
//            << multiPatchCells.objectPath() << endl;
//        multiPatchCells.write();
//
//
//        // Go through all points and remove extrusion on any cell in
//        // multiPatchCells
//        // (has to be done in separate loop since having one point on
//        // multipatches has to reset extrusion on all points of cell)
//
//        forAll(pointFaces, patchPointI)
//        {
//            if (extrudeStatus[patchPointI] != NOEXTRUDE)
//            {
//                const labelList& pFaces = pointFaces[patchPointI];
//
//                forAll(pFaces, i)
//                {
//                    label cellI =
//                        mesh.faceOwner()[pp.addressing()[pFaces[i]]];
//
//                    if (multiPatchCells.found(cellI))
//                    {
//                        if
//                        (
//                            unmarkExtrusion
//                            (
//                                patchPointI,
//                                patchDisp,
//                                patchNLayers,
//                                extrudeStatus
//                            )
//                        )
//                        {
//                            nChanged++;
//                        }
//                    }
//                }
//            }
//        }
//
//        reduce(nChanged, sumOp<label>());
//    }
//
//    Info<< "Prevented extrusion on " << nChanged
//        << " points due to multiple patch faces." << nl << endl;
//}


void Foam::autoLayerDriver::setNumLayers
(
    const labelList& patchToNLayers,
    const labelList& patchIDs,
    const indirectPrimitivePatch& pp,
    pointField& patchDisp,
    labelList& patchNLayers,
    List<extrudeMode>& extrudeStatus,
    label& nAddedCells
) const
{
    const fvMesh& mesh = meshRefiner_.mesh();

    Info<< nl << "Handling points with inconsistent layer specification ..."
        << endl;

    // Get for every point (really only necessary on patch external points)
    // the max and min of any patch faces using it.
    labelList maxLayers(patchNLayers.size(), labelMin);
    labelList minLayers(patchNLayers.size(), labelMax);

    forAll(patchIDs, i)
    {
        label patchI = patchIDs[i];

        const labelList& meshPoints = mesh.boundaryMesh()[patchI].meshPoints();

        label wantedLayers = patchToNLayers[patchI];

        forAll(meshPoints, patchPointI)
        {
            label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];

            maxLayers[ppPointI] = max(wantedLayers, maxLayers[ppPointI]);
            minLayers[ppPointI] = min(wantedLayers, minLayers[ppPointI]);
        }
    }

    syncTools::syncPointList
    (
        mesh,
        pp.meshPoints(),
        maxLayers,
        maxEqOp<label>(),
        labelMin            // null value
    );
    syncTools::syncPointList
    (
        mesh,
        pp.meshPoints(),
        minLayers,
        minEqOp<label>(),
        labelMax            // null value
    );

    // Unmark any point with different min and max
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    //label nConflicts = 0;

    forAll(maxLayers, i)
    {
        if (maxLayers[i] == labelMin || minLayers[i] == labelMax)
        {
            FatalErrorIn("setNumLayers(..)")
                << "Patchpoint:" << i << " coord:" << pp.localPoints()[i]
                << " maxLayers:" << maxLayers
                << " minLayers:" << minLayers
                << abort(FatalError);
        }
        else if (maxLayers[i] == minLayers[i])
        {
            // Ok setting.
            patchNLayers[i] = maxLayers[i];
        }
        else
        {
            // Inconsistent num layers between patch faces using point
            //if
            //(
            //    unmarkExtrusion
            //    (
            //        i,
            //        patchDisp,
            //        patchNLayers,
            //        extrudeStatus
            //    )
            //)
            //{
            //    nConflicts++;
            //}
            patchNLayers[i] = maxLayers[i];
        }
    }


    // Calculate number of cells to create
    nAddedCells = 0;
    forAll(pp.localFaces(), faceI)
    {
        const face& f = pp.localFaces()[faceI];

        // Get max of extrusion per point
        label nCells = 0;
        forAll(f, fp)
        {
            nCells = max(nCells, patchNLayers[f[fp]]);
        }

        nAddedCells += nCells;
    }
    reduce(nAddedCells, sumOp<label>());

    //reduce(nConflicts, sumOp<label>());
    //
    //Info<< "Set displacement to zero for " << nConflicts
    //    << " points due to points being on multiple regions"
    //    << " with inconsistent nLayers specification." << endl;
}


// Construct pointVectorField with correct boundary conditions for adding
// layers
Foam::tmp<Foam::pointVectorField>
Foam::autoLayerDriver::makeLayerDisplacementField
(
    const pointMesh& pMesh,
    const labelList& numLayers
)
{
    // Construct displacement field.
    const pointBoundaryMesh& pointPatches = pMesh.boundary();

    wordList patchFieldTypes
    (
        pointPatches.size(),
        slipPointPatchVectorField::typeName
    );
    wordList actualPatchTypes(patchFieldTypes.size());
    forAll(pointPatches, patchI)
    {
        actualPatchTypes[patchI] = pointPatches[patchI].type();
    }

    forAll(numLayers, patchI)
    {
        //  0 layers: do not allow slip so fixedValue 0
        // >0 layers: fixedValue which gets adapted
        if (numLayers[patchI] == 0)
        {
            patchFieldTypes[patchI] =
                zeroFixedValuePointPatchVectorField::typeName;
        }
        else if (numLayers[patchI] > 0)
        {
            patchFieldTypes[patchI] = fixedValuePointPatchVectorField::typeName;
        }
    }

    forAll(pointPatches, patchI)
    {
        if (isA<processorPointPatch>(pointPatches[patchI]))
        {
            patchFieldTypes[patchI] = calculatedPointPatchVectorField::typeName;
        }
        else if (isA<cyclicPointPatch>(pointPatches[patchI]))
        {
            patchFieldTypes[patchI] = cyclicSlipPointPatchVectorField::typeName;
        }
    }


    const polyMesh& mesh = pMesh();

    // Note: time().timeName() instead of meshRefinement::timeName() since
    // postprocessable field.
    tmp<pointVectorField> tfld
    (
        new pointVectorField
        (
            IOobject
            (
                "pointDisplacement",
                mesh.time().timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            pMesh,
            dimensionedVector("displacement", dimLength, vector::zero),
            patchFieldTypes,
            actualPatchTypes
        )
    );
    return tfld;
}


void Foam::autoLayerDriver::growNoExtrusion
(
    const indirectPrimitivePatch& pp,
    pointField& patchDisp,
    labelList& patchNLayers,
    List<extrudeMode>& extrudeStatus
) const
{
    Info<< nl << "Growing non-extrusion points by one layer ..." << endl;

    List<extrudeMode> grownExtrudeStatus(extrudeStatus);

    const faceList& localFaces = pp.localFaces();

    label nGrown = 0;

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

        bool hasSqueeze = false;
        forAll(f, fp)
        {
            if (extrudeStatus[f[fp]] == NOEXTRUDE)
            {
                hasSqueeze = true;
                break;
            }
        }

        if (hasSqueeze)
        {
            // Squeeze all points of face
            forAll(f, fp)
            {
                if
                (
                    extrudeStatus[f[fp]] == EXTRUDE
                 && grownExtrudeStatus[f[fp]] != NOEXTRUDE
                )
                {
                    grownExtrudeStatus[f[fp]] = NOEXTRUDE;
                    nGrown++;
                }
            }
        }
    }

    extrudeStatus.transfer(grownExtrudeStatus);


    // Synchronise since might get called multiple times.
    // Use the fact that NOEXTRUDE is the minimum value.
    {
        labelList status(extrudeStatus.size());
        forAll(status, i)
        {
            status[i] = extrudeStatus[i];
        }
        syncTools::syncPointList
        (
            meshRefiner_.mesh(),
            pp.meshPoints(),
            status,
            minEqOp<label>(),
            labelMax            // null value
        );
        forAll(status, i)
        {
            extrudeStatus[i] = extrudeMode(status[i]);
        }
    }


    forAll(extrudeStatus, patchPointI)
    {
        if (extrudeStatus[patchPointI] == NOEXTRUDE)
        {
            patchDisp[patchPointI] = vector::zero;
            patchNLayers[patchPointI] = 0;
        }
    }

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

    Info<< "Set displacement to zero for an additional " << nGrown
        << " points." << endl;
}


void Foam::autoLayerDriver::determineSidePatches
(
    const globalIndex& globalFaces,
    const labelListList& edgeGlobalFaces,
    const indirectPrimitivePatch& pp,

    labelList& edgePatchID,
    labelList& edgeZoneID,
    boolList& edgeFlip,
    labelList& inflateFaceID
)
{
    // Sometimes edges-to-be-extruded are on more than 2 processors.
    // Work out which 2 hold the faces to be extruded and thus which procpatch
    // the edge-face should be in. As an additional complication this might
    // mean that 2 procesors that were only edge-connected now suddenly need
    // to become face-connected i.e. have a processor patch between them.

    fvMesh& mesh = meshRefiner_.mesh();

    // Determine edgePatchID. Any additional processor boundary gets added to
    // patchToNbrProc,nbrProcToPatch and nPatches gets set to the new number
    // of patches.
    label nPatches;
    Map<label> nbrProcToPatch;
    Map<label> patchToNbrProc;
    addPatchCellLayer::calcExtrudeInfo
    (
        true,           // zoneFromAnyFace

        mesh,
        globalFaces,
        edgeGlobalFaces,
        pp,

        edgePatchID,
        nPatches,
        nbrProcToPatch,
        patchToNbrProc,
        edgeZoneID,
        edgeFlip,
        inflateFaceID
    );

    label nOldPatches = mesh.boundaryMesh().size();
    label nAdded = returnReduce(nPatches-nOldPatches, sumOp<label>());
    Info<< nl << "Adding in total " << nAdded/2 << " inter-processor patches to"
        << " handle extrusion of non-manifold processor boundaries."
        << endl;

    if (nAdded > 0)
    {
        // We might not add patches in same order as in patchToNbrProc
        // so prepare to renumber edgePatchID
        Map<label> wantedToAddedPatch;

        for (label patchI = nOldPatches; patchI < nPatches; patchI++)
        {
            label nbrProcI = patchToNbrProc[patchI];
            word name =
                    "procBoundary"
                  + Foam::name(Pstream::myProcNo())
                  + "to"
                  + Foam::name(nbrProcI);

            dictionary patchDict;
            patchDict.add("type", processorPolyPatch::typeName);
            patchDict.add("myProcNo", Pstream::myProcNo());
            patchDict.add("neighbProcNo", nbrProcI);
            patchDict.add("nFaces", 0);
            patchDict.add("startFace", mesh.nFaces());

            //Pout<< "Adding patch " << patchI
            //    << " name:" << name
            //    << " between " << Pstream::myProcNo()
            //    << " and " << nbrProcI << endl;

            label procPatchI = meshRefiner_.appendPatch
            (
                mesh,
                mesh.boundaryMesh().size(), // new patch index
                name,
                patchDict
            );
            wantedToAddedPatch.insert(patchI, procPatchI);
        }

        // Renumber edgePatchID
        forAll(edgePatchID, i)
        {
            label patchI = edgePatchID[i];
            Map<label>::const_iterator fnd = wantedToAddedPatch.find(patchI);
            if (fnd != wantedToAddedPatch.end())
            {
                edgePatchID[i] = fnd();
            }
        }

        mesh.clearOut();
        const_cast<polyBoundaryMesh&>(mesh.boundaryMesh()).updateMesh();
    }
}


void Foam::autoLayerDriver::calculateLayerThickness
(
    const indirectPrimitivePatch& pp,
    const labelList& patchIDs,
    const layerParameters& layerParams,
    const labelList& cellLevel,
    const labelList& patchNLayers,
    const scalar edge0Len,

    scalarField& thickness,
    scalarField& minThickness,
    scalarField& expansionRatio
) const
{
    const fvMesh& mesh = meshRefiner_.mesh();
    const polyBoundaryMesh& patches = mesh.boundaryMesh();


    // Rework patch-wise layer parameters into minimum per point
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Note: only layer parameters consistent with layer specification
    // method (see layerParameters) will be correct.
    scalarField firstLayerThickness(pp.nPoints(), GREAT);
    scalarField finalLayerThickness(pp.nPoints(), GREAT);
    scalarField totalThickness(pp.nPoints(), GREAT);
    scalarField expRatio(pp.nPoints(), GREAT);

    minThickness.setSize(pp.nPoints());
    minThickness = GREAT;

    forAll(patchIDs, i)
    {
        label patchI = patchIDs[i];

        const labelList& meshPoints = patches[patchI].meshPoints();

        forAll(meshPoints, patchPointI)
        {
            label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];

            firstLayerThickness[ppPointI] = min
            (
                firstLayerThickness[ppPointI],
                layerParams.firstLayerThickness()[patchI]
            );
            finalLayerThickness[ppPointI] = min
            (
                finalLayerThickness[ppPointI],
                layerParams.finalLayerThickness()[patchI]
            );
            totalThickness[ppPointI] = min
            (
                totalThickness[ppPointI],
                layerParams.thickness()[patchI]
            );
            expRatio[ppPointI] = min
            (
                expRatio[ppPointI],
                layerParams.expansionRatio()[patchI]
            );
            minThickness[ppPointI] = min
            (
                minThickness[ppPointI],
                layerParams.minThickness()[patchI]
            );
        }
    }

    syncTools::syncPointList
    (
        mesh,
        pp.meshPoints(),
        firstLayerThickness,
        minEqOp<scalar>(),
        GREAT               // null value
    );
    syncTools::syncPointList
    (
        mesh,
        pp.meshPoints(),
        finalLayerThickness,
        minEqOp<scalar>(),
        GREAT               // null value
    );
    syncTools::syncPointList
    (
        mesh,
        pp.meshPoints(),
        totalThickness,
        minEqOp<scalar>(),
        GREAT               // null value
    );
    syncTools::syncPointList
    (
        mesh,
        pp.meshPoints(),
        expRatio,
        minEqOp<scalar>(),
        GREAT               // null value
    );
    syncTools::syncPointList
    (
        mesh,
        pp.meshPoints(),
        minThickness,
        minEqOp<scalar>(),
        GREAT               // null value
    );


    // Now the thicknesses are set according to the minimum of connected
    // patches.


    // Rework relative thickness into absolute
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // by multiplying with the internal cell size.

    if (layerParams.relativeSizes())
    {
        if
        (
            min(layerParams.minThickness()) < 0
         || max(layerParams.minThickness()) > 2
        )
        {
            FatalErrorIn("calculateLayerThickness(..)")
                << "Thickness should be factor of local undistorted cell size."
                << " Valid values are [0..2]." << nl
                << " minThickness:" << layerParams.minThickness()
                << exit(FatalError);
        }


        // Determine per point the max cell level of connected cells
        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

        labelList maxPointLevel(pp.nPoints(), labelMin);

        forAll(pp, i)
        {
            label ownLevel = cellLevel[mesh.faceOwner()[pp.addressing()[i]]];

            const face& f = pp.localFaces()[i];

            forAll(f, fp)
            {
                maxPointLevel[f[fp]] = max(maxPointLevel[f[fp]], ownLevel);
            }
        }

        syncTools::syncPointList
        (
            mesh,
            pp.meshPoints(),
            maxPointLevel,
            maxEqOp<label>(),
            labelMin            // null value
        );


        forAll(maxPointLevel, pointI)
        {
            // Find undistorted edge size for this level.
            scalar edgeLen = edge0Len/(1<<maxPointLevel[pointI]);
            firstLayerThickness[pointI] *= edgeLen;
            finalLayerThickness[pointI] *= edgeLen;
            totalThickness[pointI] *= edgeLen;
            minThickness[pointI] *= edgeLen;
        }
    }



    // Rework thickness parameters into overall thickness
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    forAll(firstLayerThickness, pointI)
    {
        thickness[pointI] = layerParams.layerThickness
        (
            patchNLayers[pointI],
            firstLayerThickness[pointI],
            finalLayerThickness[pointI],
            totalThickness[pointI],
            expRatio[pointI]
        );

        expansionRatio[pointI] = layerParams.layerExpansionRatio
        (
            patchNLayers[pointI],
            firstLayerThickness[pointI],
            finalLayerThickness[pointI],
            totalThickness[pointI],
            expRatio[pointI]
        );
    }

    //Info<< "calculateLayerThickness : min:" << gMin(thickness)
    //    << " max:" << gMax(thickness) << endl;

    // Print a bit
    {
        const polyBoundaryMesh& patches = mesh.boundaryMesh();

        int oldPrecision = Info().precision();

        // Find maximum length of a patch name, for a nicer output
        label maxPatchNameLen = 0;
        forAll(patchIDs, i)
        {
            label patchI = patchIDs[i];
            word patchName = patches[patchI].name();
            maxPatchNameLen = max(maxPatchNameLen, label(patchName.size()));
        }

        Info<< nl
            << setf(ios_base::left) << setw(maxPatchNameLen) << "patch"
            << setw(0) << " faces    layers avg thickness[m]" << nl
            << setf(ios_base::left) << setw(maxPatchNameLen) << " "
            << setw(0) << "                 near-wall overall" << nl
            << setf(ios_base::left) << setw(maxPatchNameLen) << "-----"
            << setw(0) << " -----    ------ --------- -------" << endl;


        const PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh));

        forAll(patchIDs, i)
        {
            label patchI = patchIDs[i];

            const labelList& meshPoints = patches[patchI].meshPoints();

            scalar sumThickness = 0;
            scalar sumNearWallThickness = 0;
            label nMasterPoints = 0;

            forAll(meshPoints, patchPointI)
            {
                label meshPointI = meshPoints[patchPointI];
                if (isMasterPoint[meshPointI])
                {
                    label ppPointI = pp.meshPointMap()[meshPointI];

                    sumThickness += thickness[ppPointI];
                    sumNearWallThickness += layerParams.firstLayerThickness
                    (
                        patchNLayers[ppPointI],
                        firstLayerThickness[ppPointI],
                        finalLayerThickness[ppPointI],
                        thickness[ppPointI],
                        expansionRatio[ppPointI]
                    );
                    nMasterPoints++;
                }
            }

            label totNPoints = returnReduce(nMasterPoints, sumOp<label>());

            // For empty patches, totNPoints is 0.
            scalar avgThickness = 0;
            scalar avgNearWallThickness = 0;

            if (totNPoints > 0)
            {
                avgThickness =
                    returnReduce(sumThickness, sumOp<scalar>())
                  / totNPoints;
                avgNearWallThickness =
                    returnReduce(sumNearWallThickness, sumOp<scalar>())
                  / totNPoints;
            }

            Info<< setf(ios_base::left) << setw(maxPatchNameLen)
                << patches[patchI].name() << setprecision(3)
                << " " << setw(8)
                << returnReduce(patches[patchI].size(), sumOp<scalar>())
                << " " << setw(6) << layerParams.numLayers()[patchI]
                << " " << setw(8) << avgNearWallThickness
                << "  " << setw(8) << avgThickness
                << endl;
        }
        Info<< setprecision(oldPrecision) << endl;
    }
}


// Synchronize displacement among coupled patches.
void Foam::autoLayerDriver::syncPatchDisplacement
(
    const indirectPrimitivePatch& pp,
    const scalarField& minThickness,
    pointField& patchDisp,
    labelList& patchNLayers,
    List<extrudeMode>& extrudeStatus
) const
{
    const fvMesh& mesh = meshRefiner_.mesh();
    const labelList& meshPoints = pp.meshPoints();

    label nChangedTotal = 0;

    while (true)
    {
        label nChanged = 0;

        // Sync displacement (by taking min)
        syncTools::syncPointList
        (
            mesh,
            meshPoints,
            patchDisp,
            minMagSqrEqOp<vector>(),
            point::rootMax      // null value
        );

        // Unmark if displacement too small
        forAll(patchDisp, i)
        {
            if (mag(patchDisp[i]) < minThickness[i])
            {
                if
                (
                    unmarkExtrusion
                    (
                        i,
                        patchDisp,
                        patchNLayers,
                        extrudeStatus
                    )
                )
                {
                    nChanged++;
                }
            }
        }

        labelList syncPatchNLayers(patchNLayers);

        syncTools::syncPointList
        (
            mesh,
            meshPoints,
            syncPatchNLayers,
            minEqOp<label>(),
            labelMax            // null value
        );

        // Reset if differs
        // 1. take max
        forAll(syncPatchNLayers, i)
        {
            if (syncPatchNLayers[i] != patchNLayers[i])
            {
                if
                (
                    unmarkExtrusion
                    (
                        i,
                        patchDisp,
                        patchNLayers,
                        extrudeStatus
                    )
                )
                {
                    nChanged++;
                }
            }
        }

        syncTools::syncPointList
        (
            mesh,
            meshPoints,
            syncPatchNLayers,
            maxEqOp<label>(),
            labelMin            // null value
        );

        // Reset if differs
        // 2. take min
        forAll(syncPatchNLayers, i)
        {
            if (syncPatchNLayers[i] != patchNLayers[i])
            {
                if
                (
                    unmarkExtrusion
                    (
                        i,
                        patchDisp,
                        patchNLayers,
                        extrudeStatus
                    )
                )
                {
                    nChanged++;
                }
            }
        }
        nChangedTotal += nChanged;

        if (!returnReduce(nChanged, sumOp<label>()))
        {
            break;
        }
    }

    //Info<< "Prevented extrusion on "
    //    << returnReduce(nChangedTotal, sumOp<label>())
    //    << " coupled patch points during syncPatchDisplacement." << endl;
}


// Calculate displacement vector for all patch points. Uses pointNormal.
// Checks that displaced patch point would be visible from all centres
// of the faces using it.
// extrudeStatus is both input and output and gives the status of each
// patch point.
void Foam::autoLayerDriver::getPatchDisplacement
(
    const indirectPrimitivePatch& pp,
    const scalarField& thickness,
    const scalarField& minThickness,
    pointField& patchDisp,
    labelList& patchNLayers,
    List<extrudeMode>& extrudeStatus
) const
{
    Info<< nl << "Determining displacement for added points"
        << " according to pointNormal ..." << endl;

    const fvMesh& mesh = meshRefiner_.mesh();
    const vectorField& faceNormals = pp.faceNormals();
    const labelListList& pointFaces = pp.pointFaces();
    const pointField& localPoints = pp.localPoints();

    // Determine pointNormal
    // ~~~~~~~~~~~~~~~~~~~~~

    pointField pointNormals(PatchTools::pointNormals(mesh, pp));


    // Determine local length scale on patch
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    // Start off from same thickness everywhere (except where no extrusion)
    patchDisp = thickness*pointNormals;


    label nNoVisNormal = 0;
    label nExtrudeRemove = 0;


    // Check if no extrude possible.
    forAll(pointNormals, patchPointI)
    {
        label meshPointI = pp.meshPoints()[patchPointI];

        if (extrudeStatus[patchPointI] == NOEXTRUDE)
        {
            // Do not use unmarkExtrusion; forcibly set to zero extrusion.
            patchNLayers[patchPointI] = 0;
            patchDisp[patchPointI] = vector::zero;
        }
        else
        {
            // Get normal
            const vector& n = pointNormals[patchPointI];

            if (!meshTools::visNormal(n, faceNormals, pointFaces[patchPointI]))
            {
                if (debug&meshRefinement::ATTRACTION)
                {
                    Pout<< "No valid normal for point " << meshPointI
                        << ' ' << pp.points()[meshPointI]
                        << "; setting displacement to "
                        << patchDisp[patchPointI]
                        << endl;
                }

                extrudeStatus[patchPointI] = EXTRUDEREMOVE;
                nNoVisNormal++;
            }
        }
    }

    // At illegal points make displacement average of new neighbour positions
    forAll(extrudeStatus, patchPointI)
    {
        if (extrudeStatus[patchPointI] == EXTRUDEREMOVE)
        {
            point avg(vector::zero);
            label nPoints = 0;

            const labelList& pEdges = pp.pointEdges()[patchPointI];

            forAll(pEdges, i)
            {
                label edgeI = pEdges[i];

                label otherPointI = pp.edges()[edgeI].otherVertex(patchPointI);

                if (extrudeStatus[otherPointI] != NOEXTRUDE)
                {
                    avg += localPoints[otherPointI] + patchDisp[otherPointI];
                    nPoints++;
                }
            }

            if (nPoints > 0)
            {
                if (debug&meshRefinement::ATTRACTION)
                {
                    Pout<< "Displacement at illegal point "
                        << localPoints[patchPointI]
                        << " set to "
                        << (avg / nPoints - localPoints[patchPointI])
                        << endl;
                }

                patchDisp[patchPointI] =
                    avg / nPoints
                  - localPoints[patchPointI];

                nExtrudeRemove++;
            }
            else
            {
                // All surrounding points are not extruded. Leave patchDisp
                // intact.
            }
        }
    }

    Info<< "Detected " << returnReduce(nNoVisNormal, sumOp<label>())
        << " points with point normal pointing through faces." << nl
        << "Reset displacement at "
        << returnReduce(nExtrudeRemove, sumOp<label>())
        << " points to average of surrounding points." << endl;

    // Make sure displacement is equal on both sides of coupled patches.
    syncPatchDisplacement
    (
        pp,
        minThickness,
        patchDisp,
        patchNLayers,
        extrudeStatus
    );

    Info<< endl;
}


bool Foam::autoLayerDriver::sameEdgeNeighbour
(
    const labelListList& globalEdgeFaces,
    const label myGlobalFaceI,
    const label nbrGlobFaceI,
    const label edgeI
) const
{
    const labelList& eFaces = globalEdgeFaces[edgeI];
    if (eFaces.size() == 2)
    {
        return edge(myGlobalFaceI, nbrGlobFaceI) == edge(eFaces[0], eFaces[1]);
    }
    else
    {
        return false;
    }
}


void Foam::autoLayerDriver::getVertexString
(
    const indirectPrimitivePatch& pp,
    const labelListList& globalEdgeFaces,
    const label faceI,
    const label edgeI,
    const label myGlobFaceI,
    const label nbrGlobFaceI,
    DynamicList<label>& vertices
) const
{
    const labelList& fEdges = pp.faceEdges()[faceI];
    label fp = findIndex(fEdges, edgeI);

    if (fp == -1)
    {
        FatalErrorIn("autoLayerDriver::getVertexString(..)")
            << "problem." << abort(FatalError);
    }

    // Search back
    label startFp = fp;

    forAll(fEdges, i)
    {
        label prevFp = fEdges.rcIndex(startFp);
        if
        (
           !sameEdgeNeighbour
            (
                globalEdgeFaces,
                myGlobFaceI,
                nbrGlobFaceI,
                fEdges[prevFp]
            )
        )
        {
            break;
        }
        startFp = prevFp;
    }

    label endFp = fp;
    forAll(fEdges, i)
    {
        label nextFp = fEdges.fcIndex(endFp);
        if
        (
           !sameEdgeNeighbour
            (
                globalEdgeFaces,
                myGlobFaceI,
                nbrGlobFaceI,
                fEdges[nextFp]
            )
        )
        {
            break;
        }
        endFp = nextFp;
    }

    const face& f = pp.localFaces()[faceI];
    vertices.clear();
    fp = startFp;
    while (fp != endFp)
    {
        vertices.append(f[fp]);
        fp = f.fcIndex(fp);
    }
    vertices.append(f[fp]);
    fp = f.fcIndex(fp);
    vertices.append(f[fp]);
}


// Truncates displacement
// - for all patchFaces in the faceset displacement gets set to zero
// - all displacement < minThickness gets set to zero
Foam::label Foam::autoLayerDriver::truncateDisplacement
(
    const globalIndex& globalFaces,
    const labelListList& edgeGlobalFaces,
    const indirectPrimitivePatch& pp,
    const scalarField& minThickness,
    const faceSet& illegalPatchFaces,
    pointField& patchDisp,
    labelList& patchNLayers,
    List<extrudeMode>& extrudeStatus
) const
{
    const fvMesh& mesh = meshRefiner_.mesh();

    label nChanged = 0;

    const Map<label>& meshPointMap = pp.meshPointMap();

    forAllConstIter(faceSet, illegalPatchFaces, iter)
    {
        label faceI = iter.key();

        if (mesh.isInternalFace(faceI))
        {
            FatalErrorIn("truncateDisplacement(..)")
                << "Faceset " << illegalPatchFaces.name()
                << " contains internal face " << faceI << nl
                << "It should only contain patch faces" << abort(FatalError);
        }

        const face& f = mesh.faces()[faceI];


        forAll(f, fp)
        {
            if (meshPointMap.found(f[fp]))
            {
                label patchPointI = meshPointMap[f[fp]];

                if (extrudeStatus[patchPointI] != NOEXTRUDE)
                {
                    unmarkExtrusion
                    (
                        patchPointI,
                        patchDisp,
                        patchNLayers,
                        extrudeStatus
                    );
                    nChanged++;
                }
            }
        }
    }

    forAll(patchDisp, patchPointI)
    {
        if (mag(patchDisp[patchPointI]) < minThickness[patchPointI])
        {
            if
            (
                unmarkExtrusion
                (
                    patchPointI,
                    patchDisp,
                    patchNLayers,
                    extrudeStatus
                )
            )
            {
                nChanged++;
            }
        }
        else if (extrudeStatus[patchPointI] == NOEXTRUDE)
        {
            // Make sure displacement is 0. Should already be so but ...
            patchDisp[patchPointI] = vector::zero;
            patchNLayers[patchPointI] = 0;
        }
    }


    const faceList& localFaces = pp.localFaces();

    while (true)
    {
        syncPatchDisplacement
        (
            pp,
            minThickness,
            patchDisp,
            patchNLayers,
            extrudeStatus
        );


        // Pinch
        // ~~~~~

        // Make sure that a face doesn't have two non-consecutive areas
        // not extruded (e.g. quad where vertex 0 and 2 are not extruded
        // but 1 and 3 are) since this gives topological errors.

        label nPinched = 0;

        forAll(localFaces, i)
        {
            const face& localF = localFaces[i];

            // Count number of transitions from unsnapped to snapped.
            label nTrans = 0;

            extrudeMode prevMode = extrudeStatus[localF.prevLabel(0)];

            forAll(localF, fp)
            {
                extrudeMode fpMode = extrudeStatus[localF[fp]];

                if (prevMode == NOEXTRUDE && fpMode != NOEXTRUDE)
                {
                    nTrans++;
                }
                prevMode = fpMode;
            }

            if (nTrans > 1)
            {
                // Multiple pinches. Reset whole face as unextruded.
                if
                (
                    unmarkExtrusion
                    (
                        localF,
                        patchDisp,
                        patchNLayers,
                        extrudeStatus
                    )
                )
                {
                    nPinched++;
                    nChanged++;
                }
            }
        }

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

        Info<< "truncateDisplacement : Unextruded " << nPinched
            << " faces due to non-consecutive vertices being extruded." << endl;


        // Butterfly
        // ~~~~~~~~~

        // Make sure that a string of edges becomes a single face so
        // not a butterfly. Occassionally an 'edge' will have a single dangling
        // vertex due to face combining. These get extruded as a single face
        // (with a dangling vertex) so make sure this extrusion forms a single
        // shape.
        //  - continuous i.e. no butterfly:
        //      +     +
        //      |\   /|
        //      | \ / |
        //      +--+--+
        //  - extrudes from all but the endpoints i.e. no partial
        //    extrude
        //            +
        //           /|
        //          / |
        //      +--+--+
        // The common error topology is a pinch somewhere in the middle
        label nButterFly = 0;
        {
            DynamicList<label> stringedVerts;
            forAll(pp.edges(), edgeI)
            {
                const labelList& globFaces = edgeGlobalFaces[edgeI];

                if (globFaces.size() == 2)
                {
                    label myFaceI = pp.edgeFaces()[edgeI][0];
                    label myGlobalFaceI = globalFaces.toGlobal
                    (
                        pp.addressing()[myFaceI]
                    );
                    label nbrGlobalFaceI =
                    (
                        globFaces[0] != myGlobalFaceI
                      ? globFaces[0]
                      : globFaces[1]
                    );
                    getVertexString
                    (
                        pp,
                        edgeGlobalFaces,
                        myFaceI,
                        edgeI,
                        myGlobalFaceI,
                        nbrGlobalFaceI,
                        stringedVerts
                    );

                    if
                    (
                        extrudeStatus[stringedVerts[0]] != NOEXTRUDE
                     || extrudeStatus[stringedVerts.last()] != NOEXTRUDE
                    )
                    {
                        // Any pinch in the middle
                        bool pinch = false;
                        for (label i = 1; i < stringedVerts.size()-1; i++)
                        {
                            if
                            (
                                extrudeStatus[stringedVerts[i]] == NOEXTRUDE
                            )
                            {
                                pinch = true;
                                break;
                            }
                        }
                        if (pinch)
                        {
                            forAll(stringedVerts, i)
                            {
                                if
                                (
                                    unmarkExtrusion
                                    (
                                        stringedVerts[i],
                                        patchDisp,
                                        patchNLa