autoSnapDriver.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2015 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 snapping to the surface
 
\*----------------------------------------------------------------------------*/
 
#include "autoSnapDriver.H"
#include "motionSmoother.H"
#include "polyTopoChange.H"
#include "syncTools.H"
#include "fvMesh.H"
#include "Time.H"
#include "OFstream.H"
#include "OBJstream.H"
#include "mapPolyMesh.H"
#include "pointEdgePoint.H"
#include "PointEdgeWave.H"
#include "mergePoints.H"
#include "snapParameters.H"
#include "refinementSurfaces.H"
#include "searchableSurfaces.H"
#include "unitConversion.H"
#include "localPointRegion.H"
#include "PatchTools.H"
#include "refinementFeatures.H"
 
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
namespace Foam
{
 
defineTypeNameAndDebug(autoSnapDriver, 0);
 
} // End namespace Foam
 
 
// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
// Calculate geometrically collocated points, Requires PackedList to be
// sized and initalised!
Foam::label Foam::autoSnapDriver::getCollocatedPoints
(
    const scalar tol,
    const pointField& points,
    PackedBoolList& isCollocatedPoint
)
{
    labelList pointMap;
    label nUnique = mergePoints
    (
        points,                         // points
        tol,                            // mergeTol
        false,                          // verbose
        pointMap
    );
    bool hasMerged = (nUnique < points.size());
 
    if (!returnReduce(hasMerged, orOp<bool>()))
    {
        return 0;
    }
 
    // Determine which merged points are referenced more than once
    label nCollocated = 0;
 
    // Per old point the newPoint. Or -1 (not set yet) or -2 (already seen
    // twice)
    labelList firstOldPoint(nUnique, -1);
    forAll(pointMap, oldPointI)
    {
        label newPointI = pointMap[oldPointI];
 
        if (firstOldPoint[newPointI] == -1)
        {
            // First use of oldPointI. Store.
            firstOldPoint[newPointI] = oldPointI;
        }
        else if (firstOldPoint[newPointI] == -2)
        {
            // Third or more reference of oldPointI -> non-manifold
            isCollocatedPoint.set(oldPointI, 1u);
            nCollocated++;
        }
        else
        {
            // Second reference of oldPointI -> non-manifold
            isCollocatedPoint.set(firstOldPoint[newPointI], 1u);
            nCollocated++;
 
            isCollocatedPoint.set(oldPointI, 1u);
            nCollocated++;
 
            // Mark with special value to save checking next time round
            firstOldPoint[newPointI] = -2;
        }
    }
    return returnReduce(nCollocated, sumOp<label>());
}
 
 
Foam::tmp<Foam::pointField> Foam::autoSnapDriver::smoothInternalDisplacement
(
    const meshRefinement& meshRefiner,
    const motionSmoother& meshMover
)
{
    const indirectPrimitivePatch& pp = meshMover.patch();
    const polyMesh& mesh = meshMover.mesh();
 
    // Get neighbour refinement
    const hexRef8& cutter = meshRefiner.meshCutter();
    const labelList& cellLevel = cutter.cellLevel();
 
 
    // Get the faces on the boundary
    PackedBoolList isFront(mesh.nFaces());
    forAll(pp.addressing(), i)
    {
        isFront[pp.addressing()[i]] = true;
    }
 
    // Walk out from the surface a bit. Poor man's FaceCellWave.
    // Commented out for now - not sure if needed and if so how much
    //for (label iter = 0; iter < 2; iter++)
    //{
    //    PackedBoolList newIsFront(mesh.nFaces());
    //
    //    forAll(isFront, faceI)
    //    {
    //        if (isFront[faceI])
    //        {
    //            label own = mesh.faceOwner()[faceI];
    //            const cell& ownFaces = mesh.cells()[own];
    //            forAll(ownFaces, i)
    //            {
    //                newIsFront[ownFaces[i]] = true;
    //            }
    //
    //            if (mesh.isInternalFace(faceI))
    //            {
    //                label nei = mesh.faceNeighbour()[faceI];
    //                const cell& neiFaces = mesh.cells()[nei];
    //                forAll(neiFaces, i)
    //                {
    //                    newIsFront[neiFaces[i]] = true;
    //                }
    //            }
    //        }
    //    }
    //
    //    syncTools::syncFaceList
    //    (
    //        mesh,
    //        newIsFront,
    //        orEqOp<unsigned int>()
    //    );
    //
    //    isFront = newIsFront;
    //}
 
    // Mark all points on faces
    //  - not on the boundary
    //  - inbetween differing refinement levels
    PackedBoolList isMovingPoint(mesh.nPoints());
 
    label nInterface = 0;
 
    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
    {
        label ownLevel = cellLevel[mesh.faceOwner()[faceI]];
        label neiLevel = cellLevel[mesh.faceNeighbour()[faceI]];
 
        if (!isFront[faceI] && ownLevel != neiLevel)
        {
            const face& f = mesh.faces()[faceI];
            forAll(f, fp)
            {
                isMovingPoint[f[fp]] = true;
            }
 
            nInterface++;
        }
    }
 
    labelList neiCellLevel;
    syncTools::swapBoundaryCellList(mesh, cellLevel, neiCellLevel);
 
    for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
    {
        label ownLevel = cellLevel[mesh.faceOwner()[faceI]];
        label neiLevel = neiCellLevel[faceI-mesh.nInternalFaces()];
 
        if (!isFront[faceI] && ownLevel != neiLevel)
        {
            const face& f = mesh.faces()[faceI];
            forAll(f, fp)
            {
                isMovingPoint[f[fp]] = true;
            }
 
            nInterface++;
        }
    }
 
    if (debug)
    {
        reduce(nInterface, sumOp<label>());
        Info<< "Found " << nInterface << " faces out of "
            << mesh.globalData().nTotalFaces()
            << " inbetween refinement regions." << endl;
    }
 
    // Make sure that points that are coupled to a moving point are marked
    // as well
    syncTools::syncPointList(mesh, isMovingPoint, maxEqOp<unsigned int>(), 0);
 
    // Unmark any point on the boundary. If we're doing zero iterations of
    // face-cell wave we might have coupled points not being unmarked.
    forAll(pp.meshPoints(), pointI)
    {
        isMovingPoint[pp.meshPoints()[pointI]] = false;
    }
 
    // Make sure that points that are coupled to meshPoints but not on a patch
    // are unmarked as well
    syncTools::syncPointList(mesh, isMovingPoint, minEqOp<unsigned int>(), 1);
 
 
    // Calculate average of connected cells
    labelList nCells(mesh.nPoints(), 0);
    pointField sumLocation(mesh.nPoints(), vector::zero);
 
    forAll(isMovingPoint, pointI)
    {
        if (isMovingPoint[pointI])
        {
            const labelList& pCells = mesh.pointCells(pointI);
 
            forAll(pCells, i)
            {
                sumLocation[pointI] += mesh.cellCentres()[pCells[i]];
                nCells[pointI]++;
            }
        }
    }
 
    // Sum
    syncTools::syncPointList(mesh, nCells, plusEqOp<label>(), label(0));
    syncTools::syncPointList
    (
        mesh,
        sumLocation,
        plusEqOp<point>(),
        vector::zero
    );
 
    tmp<pointField> tdisplacement(new pointField(mesh.nPoints(), vector::zero));
    pointField& displacement = tdisplacement();
 
    label nAdapted = 0;
 
    forAll(displacement, pointI)
    {
        if (nCells[pointI] > 0)
        {
            displacement[pointI] =
                sumLocation[pointI]/nCells[pointI]-mesh.points()[pointI];
            nAdapted++;
        }
    }
 
    reduce(nAdapted, sumOp<label>());
    Info<< "Smoothing " << nAdapted << " points inbetween refinement regions."
        << endl;
 
    return tdisplacement;
}
 
 
// Calculate displacement as average of patch points.
Foam::tmp<Foam::pointField> Foam::autoSnapDriver::smoothPatchDisplacement
(
    const motionSmoother& meshMover,
    const List<labelPair>& baffles
)
{
    const indirectPrimitivePatch& pp = meshMover.patch();
 
    // Calculate geometrically non-manifold points on the patch to be moved.
    PackedBoolList nonManifoldPoint(pp.nPoints());
    label nNonManifoldPoints = getCollocatedPoints
    (
        SMALL,
        pp.localPoints(),
        nonManifoldPoint
    );
    Info<< "Found " << nNonManifoldPoints << " non-manifold point(s)."
        << endl;
 
 
    // Average points
    // ~~~~~~~~~~~~~~
 
    // We determine three points:
    // - average of (centres of) connected patch faces
    // - average of (centres of) connected internal mesh faces
    // - as fallback: centre of any connected cell
    // so we can do something moderately sensible for non/manifold points.
 
    // Note: the averages are calculated properly parallel. This is
    // necessary to get the points shared by processors correct.
 
 
    const labelListList& pointFaces = pp.pointFaces();
    const labelList& meshPoints = pp.meshPoints();
    const pointField& points = pp.points();
    const polyMesh& mesh = meshMover.mesh();
 
    // Get labels of faces to count (master of coupled faces and baffle pairs)
    PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh));
 
    {
        forAll(baffles, i)
        {
            label f0 = baffles[i].first();
            label f1 = baffles[i].second();
 
            if (isMasterFace.get(f0))
            {
                // Make f1 a slave
                isMasterFace.unset(f1);
            }
            else if (isMasterFace.get(f1))
            {
                isMasterFace.unset(f0);
            }
            else
            {
                FatalErrorInFunction
                    << "Both sides of baffle consisting of faces " << f0
                    << " and " << f1 << " are already slave faces."
                    << abort(FatalError);
            }
        }
    }
 
 
    // Get average position of boundary face centres
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
    vectorField avgBoundary(pointFaces.size(), vector::zero);
    labelList nBoundary(pointFaces.size(), 0);
 
    forAll(pointFaces, patchPointI)
    {
        const labelList& pFaces = pointFaces[patchPointI];
 
        forAll(pFaces, pfI)
        {
            label faceI = pFaces[pfI];
 
            if (isMasterFace.get(pp.addressing()[faceI]))
            {
                avgBoundary[patchPointI] += pp[faceI].centre(points);
                nBoundary[patchPointI]++;
            }
        }
    }
 
    syncTools::syncPointList
    (
        mesh,
        pp.meshPoints(),
        avgBoundary,
        plusEqOp<point>(),  // combine op
        vector::zero        // null value
    );
    syncTools::syncPointList
    (
        mesh,
        pp.meshPoints(),
        nBoundary,
        plusEqOp<label>(),  // combine op
        label(0)            // null value
    );
 
    forAll(avgBoundary, i)
    {
        avgBoundary[i] /= nBoundary[i];
    }
 
 
    // Get average position of internal face centres
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
    vectorField avgInternal;
    labelList nInternal;
    {
        vectorField globalSum(mesh.nPoints(), vector::zero);
        labelList globalNum(mesh.nPoints(), 0);
 
        // Note: no use of pointFaces
        const faceList& faces = mesh.faces();
 
        for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
        {
            const face& f = faces[faceI];
            const point& fc = mesh.faceCentres()[faceI];
 
            forAll(f, fp)
            {
                globalSum[f[fp]] += fc;
                globalNum[f[fp]]++;
            }
        }
 
        // Count coupled faces as internal ones (but only once)
        const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
        forAll(patches, patchI)
        {
            if
            (
                patches[patchI].coupled()
             && refCast<const coupledPolyPatch>(patches[patchI]).owner()
            )
            {
                const coupledPolyPatch& pp =
                    refCast<const coupledPolyPatch>(patches[patchI]);
 
                const vectorField::subField faceCentres = pp.faceCentres();
 
                forAll(pp, i)
                {
                    const face& f = pp[i];
                    const point& fc = faceCentres[i];
 
                    forAll(f, fp)
                    {
                        globalSum[f[fp]] += fc;
                        globalNum[f[fp]]++;
                    }
                }
            }
        }
 
        syncTools::syncPointList
        (
            mesh,
            globalSum,
            plusEqOp<vector>(), // combine op
            vector::zero        // null value
        );
        syncTools::syncPointList
        (
            mesh,
            globalNum,
            plusEqOp<label>(),  // combine op
            label(0)            // null value
        );
 
        avgInternal.setSize(meshPoints.size());
        nInternal.setSize(meshPoints.size());
 
        forAll(avgInternal, patchPointI)
        {
            label meshPointI = meshPoints[patchPointI];
 
            nInternal[patchPointI] = globalNum[meshPointI];
 
            if (nInternal[patchPointI] == 0)
            {
                avgInternal[patchPointI] = globalSum[meshPointI];
            }
            else
            {
                avgInternal[patchPointI] =
                    globalSum[meshPointI]
                  / nInternal[patchPointI];
            }
        }
    }
 
 
    // Precalculate any cell using mesh point (replacement of pointCells()[])
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
    labelList anyCell(mesh.nPoints(), -1);
    forAll(mesh.faceOwner(), faceI)
    {
        label own = mesh.faceOwner()[faceI];
        const face& f = mesh.faces()[faceI];
 
        forAll(f, fp)
        {
            anyCell[f[fp]] = own;
        }
    }
 
 
    // Displacement to calculate.
    tmp<pointField> tpatchDisp(new pointField(meshPoints.size(), vector::zero));
    pointField& patchDisp = tpatchDisp();
 
    forAll(pointFaces, i)
    {
        label meshPointI = meshPoints[i];
        const point& currentPos = pp.points()[meshPointI];
 
        // Now we have the two average points: avgBoundary and avgInternal
        // and how many boundary/internal faces connect to the point
        // (nBoundary, nInternal)
        // Do some blending between the two.
        // Note: the following section has some reasoning behind it but the
        // blending factors can be experimented with.
 
        point newPos;
 
        if (!nonManifoldPoint.get(i))
        {
            // Points that are manifold. Weight the internal and boundary
            // by their number of faces and blend with
            scalar internalBlend = 0.1;
            scalar blend = 0.1;
 
            point avgPos =
                (
                   internalBlend*nInternal[i]*avgInternal[i]
                  +(1-internalBlend)*nBoundary[i]*avgBoundary[i]
                )
              / (internalBlend*nInternal[i]+(1-internalBlend)*nBoundary[i]);
 
            newPos = (1-blend)*avgPos + blend*currentPos;
        }
        else if (nInternal[i] == 0)
        {
            // Non-manifold without internal faces. Use any connected cell
            // as internal point instead. Use precalculated any cell to avoid
            // e.g. pointCells()[meshPointI][0]
 
            const point& cc = mesh.cellCentres()[anyCell[meshPointI]];
 
            scalar cellCBlend = 0.8;
            scalar blend = 0.1;
 
            point avgPos = (1-cellCBlend)*avgBoundary[i] + cellCBlend*cc;
 
            newPos = (1-blend)*avgPos + blend*currentPos;
        }
        else
        {
            // Non-manifold point with internal faces connected to them
            scalar internalBlend = 0.9;
            scalar blend = 0.1;
 
            point avgPos =
                internalBlend*avgInternal[i]
              + (1-internalBlend)*avgBoundary[i];
 
            newPos = (1-blend)*avgPos + blend*currentPos;
        }
 
        patchDisp[i] = newPos - currentPos;
    }
 
    return tpatchDisp;
}
//XXXXXXX
//Foam::tmp<Foam::pointField> Foam::autoSnapDriver::avg
//(
//    const indirectPrimitivePatch& pp,
//    const pointField& localPoints
//)
//{
//    const labelListList& pointEdges = pp.pointEdges();
//    const edgeList& edges = pp.edges();
//
//    tmp<pointField> tavg(new pointField(pointEdges.size(), vector::zero));
//    pointField& avg = tavg();
//
//    forAll(pointEdges, vertI)
//    {
//        vector& avgPos = avg[vertI];
//
//        const labelList& pEdges = pointEdges[vertI];
//
//        forAll(pEdges, myEdgeI)
//        {
//            const edge& e = edges[pEdges[myEdgeI]];
//
//            label otherVertI = e.otherVertex(vertI);
//
//            avgPos += localPoints[otherVertI];
//        }
//
//        avgPos /= pEdges.size();
//    }
//    return tavg;
//}
//Foam::tmp<Foam::pointField>
//Foam::autoSnapDriver::smoothLambdaMuPatchDisplacement
//(
//    const motionSmoother& meshMover,
//    const List<labelPair>& baffles
//)
//{
//    const indirectPrimitivePatch& pp = meshMover.patch();
//    pointField newLocalPoints(pp.localPoints());
//
//    const label iters = 90;
//    const scalar lambda = 0.33;
//    const scalar mu = 0.34;
//
//    for (label iter = 0; iter < iters; iter++)
//    {
//        // Lambda
//        newLocalPoints =
//            (1 - lambda)*newLocalPoints
//          + lambda*avg(pp, newLocalPoints);
//
//        // Mu
//        newLocalPoints =
//            (1 + mu)*newLocalPoints
//          - mu*avg(pp, newLocalPoints);
//    }
//    return newLocalPoints-pp.localPoints();
//}
//XXXXXXX
 
 
Foam::tmp<Foam::scalarField> Foam::autoSnapDriver::edgePatchDist
(
    const pointMesh& pMesh,
    const indirectPrimitivePatch& pp
)
{
    const polyMesh& mesh = pMesh();
 
    // Set initial changed points to all the patch points
    List<pointEdgePoint> wallInfo(pp.nPoints());
 
    forAll(pp.localPoints(), ppI)
    {
        wallInfo[ppI] = pointEdgePoint(pp.localPoints()[ppI], 0.0);
    }
 
    // Current info on points
    List<pointEdgePoint> allPointInfo(mesh.nPoints());
 
    // Current info on edges
    List<pointEdgePoint> allEdgeInfo(mesh.nEdges());
 
    PointEdgeWave<pointEdgePoint> wallCalc
    (
        mesh,
        pp.meshPoints(),
        wallInfo,
 
        allPointInfo,
        allEdgeInfo,
        mesh.globalData().nTotalPoints()  // max iterations
    );
 
    // Copy edge values into scalarField
    tmp<scalarField> tedgeDist(new scalarField(mesh.nEdges()));
    scalarField& edgeDist = tedgeDist();
 
    forAll(allEdgeInfo, edgeI)
    {
        edgeDist[edgeI] = Foam::sqrt(allEdgeInfo[edgeI].distSqr());
    }
 
    return tedgeDist;
}
 
 
void Foam::autoSnapDriver::dumpMove
(
    const fileName& fName,
    const pointField& meshPts,
    const pointField& surfPts
)
{
    // Dump direction of growth into file
    Info<< "Dumping move direction to " << fName << endl;
 
    OFstream nearestStream(fName);
 
    label vertI = 0;
 
    forAll(meshPts, ptI)
    {
        meshTools::writeOBJ(nearestStream, meshPts[ptI]);
        vertI++;
 
        meshTools::writeOBJ(nearestStream, surfPts[ptI]);
        vertI++;
 
        nearestStream<< "l " << vertI-1 << ' ' << vertI << nl;
    }
}
 
 
// Check whether all displacement vectors point outwards of patch. Return true
// if so.
bool Foam::autoSnapDriver::outwardsDisplacement
(
    const indirectPrimitivePatch& pp,
    const vectorField& patchDisp
)
{
    const vectorField& faceNormals = pp.faceNormals();
    const labelListList& pointFaces = pp.pointFaces();
 
    forAll(pointFaces, pointI)
    {
        const labelList& pFaces = pointFaces[pointI];
 
        vector disp(patchDisp[pointI]);
 
        scalar magDisp = mag(disp);
 
        if (magDisp > SMALL)
        {
            disp /= magDisp;
 
            bool outwards = meshTools::visNormal(disp, faceNormals, pFaces);
 
            if (!outwards)
            {
                Warning<< "Displacement " << patchDisp[pointI]
                    << " at mesh point " << pp.meshPoints()[pointI]
                    << " coord " << pp.points()[pp.meshPoints()[pointI]]
                    << " points through the surrounding patch faces" << endl;
                return false;
            }
        }
        else
        {
            //? Displacement small but in wrong direction. Would probably be ok.
        }
    }
    return true;
}
 
 
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
Foam::autoSnapDriver::autoSnapDriver
(
    meshRefinement& meshRefiner,
    const labelList& globalToMasterPatch,
    const labelList& globalToSlavePatch
)
:
    meshRefiner_(meshRefiner),
    globalToMasterPatch_(globalToMasterPatch),
    globalToSlavePatch_(globalToSlavePatch)
{}
 
 
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
Foam::scalarField Foam::autoSnapDriver::calcSnapDistance
(
    const fvMesh& mesh,
    const snapParameters& snapParams,
    const indirectPrimitivePatch& pp
)
{
    const edgeList& edges = pp.edges();
    const labelListList& pointEdges = pp.pointEdges();
    const pointField& localPoints = pp.localPoints();
 
    scalarField maxEdgeLen(localPoints.size(), -GREAT);
 
    forAll(pointEdges, pointI)
    {
        const labelList& pEdges = pointEdges[pointI];
 
        forAll(pEdges, pEdgeI)
        {
            const edge& e = edges[pEdges[pEdgeI]];
 
            scalar len = e.mag(localPoints);
 
            maxEdgeLen[pointI] = max(maxEdgeLen[pointI], len);
        }
    }
 
    syncTools::syncPointList
    (
        mesh,
        pp.meshPoints(),
        maxEdgeLen,
        maxEqOp<scalar>(),  // combine op
        -GREAT              // null value
    );
 
    return scalarField(snapParams.snapTol()*maxEdgeLen);
}
 
 
void Foam::autoSnapDriver::preSmoothPatch
(
    const meshRefinement& meshRefiner,
    const snapParameters& snapParams,
    const label nInitErrors,
    const List<labelPair>& baffles,
    motionSmoother& meshMover
)
{
    const fvMesh& mesh = meshRefiner.mesh();
 
    labelList checkFaces;
 
    if (snapParams.nSmoothInternal() > 0)
    {
        Info<< "Smoothing patch and internal points ..." << endl;
    }
    else
    {
        Info<< "Smoothing patch points ..." << endl;
    }
 
    for
    (
        label smoothIter = 0;
        smoothIter < snapParams.nSmoothPatch();
        smoothIter++
    )
    {
        Info<< "Smoothing iteration " << smoothIter << endl;
        checkFaces.setSize(mesh.nFaces());
        forAll(checkFaces, faceI)
        {
            checkFaces[faceI] = faceI;
        }
 
        // If enabled smooth the internal points
        if (snapParams.nSmoothInternal() > smoothIter)
        {
            // Override values on internal points on refinement interfaces
            meshMover.pointDisplacement().internalField() =
                smoothInternalDisplacement(meshRefiner, meshMover);
        }
 
        // Smooth the patch points
        pointField patchDisp(smoothPatchDisplacement(meshMover, baffles));
        //pointField patchDisp
        //(
        //  smoothLambdaMuPatchDisplacement(meshMover, baffles)
        //);
 
        // Take over patch displacement as boundary condition on
        // pointDisplacement
        meshMover.setDisplacement(patchDisp);
 
        // Start off from current mesh.points()
        meshMover.correct();
 
        scalar oldErrorReduction = -1;
 
        for (label snapIter = 0; snapIter < 2*snapParams.nSnap(); snapIter++)
        {
            Info<< nl << "Scaling iteration " << snapIter << endl;
 
            if (snapIter == snapParams.nSnap())
            {
                Info<< "Displacement scaling for error reduction set to 0."
                    << endl;
                oldErrorReduction = meshMover.setErrorReduction(0.0);
            }
 
            // Try to adapt mesh to obtain displacement by smoothly
            // decreasing displacement at error locations.
            if (meshMover.scaleMesh(checkFaces, baffles, true, nInitErrors))
            {
                Info<< "Successfully moved mesh" << endl;
                break;
            }
        }
 
        if (oldErrorReduction >= 0)
        {
            meshMover.setErrorReduction(oldErrorReduction);
        }
        Info<< endl;
    }
 
 
    // The current mesh is the starting mesh to smooth from.
    meshMover.correct();
 
    if (debug&meshRefinement::MESH)
    {
        const_cast<Time&>(mesh.time())++;
        Info<< "Writing patch smoothed mesh to time "
            << meshRefiner.timeName() << '.' << endl;
        meshRefiner.write
        (
            meshRefinement::debugType(debug),
            meshRefinement::writeType
            (
                meshRefinement::writeLevel()
              | meshRefinement::WRITEMESH
            ),
            mesh.time().path()/meshRefiner.timeName()
        );
        Info<< "Dumped mesh in = "
            << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
    }
 
    Info<< "Patch points smoothed in = "
        << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
}
 
 
// Get (pp-local) indices of points that are both on zone and on patched surface
Foam::labelList Foam::autoSnapDriver::getZoneSurfacePoints
(
    const fvMesh& mesh,
    const indirectPrimitivePatch& pp,
    const word& zoneName
)
{
    label zoneI = mesh.faceZones().findZoneID(zoneName);
 
    if (zoneI == -1)
    {
        FatalErrorInFunction
            << "Cannot find zone " << zoneName
            << exit(FatalError);
    }
 
    const faceZone& fZone = mesh.faceZones()[zoneI];
 
 
    // Could use PrimitivePatch & localFaces to extract points but might just
    // as well do it ourselves.
 
    boolList pointOnZone(pp.nPoints(), false);
 
    forAll(fZone, i)
    {
        const face& f = mesh.faces()[fZone[i]];
 
        forAll(f, fp)
        {
            label meshPointI = f[fp];
 
            Map<label>::const_iterator iter =
                pp.meshPointMap().find(meshPointI);
 
            if (iter != pp.meshPointMap().end())
            {
                label pointI = iter();
                pointOnZone[pointI] = true;
            }
        }
    }
 
    return findIndices(pointOnZone, true);
}
 
 
Foam::tmp<Foam::pointField> Foam::autoSnapDriver::avgCellCentres
(
    const fvMesh& mesh,
    const indirectPrimitivePatch& pp
)
{
    const labelListList& pointFaces = pp.pointFaces();
 
 
    tmp<pointField> tavgBoundary
    (
        new pointField(pointFaces.size(), vector::zero)
    );
    pointField& avgBoundary = tavgBoundary();
    labelList nBoundary(pointFaces.size(), 0);
 
    forAll(pointFaces, pointI)
    {
        const labelList& pFaces = pointFaces[pointI];
 
        forAll(pFaces, pfI)
        {
            label faceI = pFaces[pfI];
            label meshFaceI = pp.addressing()[faceI];
 
            label own = mesh.faceOwner()[meshFaceI];
            avgBoundary[pointI] += mesh.cellCentres()[own];
            nBoundary[pointI]++;
        }
    }
 
    syncTools::syncPointList
    (
        mesh,
        pp.meshPoints(),
        avgBoundary,
        plusEqOp<point>(),  // combine op
        vector::zero        // null value
    );
    syncTools::syncPointList
    (
        mesh,
        pp.meshPoints(),
        nBoundary,
        plusEqOp<label>(),  // combine op
        label(0)            // null value
    );
 
    forAll(avgBoundary, i)
    {
        avgBoundary[i] /= nBoundary[i];
    }
    return tavgBoundary;
}
 
 
//Foam::tmp<Foam::scalarField> Foam::autoSnapDriver::calcEdgeLen
//(
//    const indirectPrimitivePatch& pp
//) const
//{
//    // Get local edge length based on refinement level
//    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//    // (Ripped from autoLayerDriver)
//
//    tmp<scalarField> tedgeLen(new scalarField(pp.nPoints()));
//    scalarField& edgeLen = tedgeLen();
//    {
//        const fvMesh& mesh = meshRefiner_.mesh();
//        const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength();
//        const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
//
//        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.
//            edgeLen[pointI] = edge0Len/(1<<maxPointLevel[pointI]);
//        }
//    }
//    return tedgeLen;
//}
 
 
void Foam::autoSnapDriver::detectNearSurfaces
(
    const scalar planarCos,
    const indirectPrimitivePatch& pp,
    const pointField& nearestPoint,
    const vectorField& nearestNormal,
 
    vectorField& disp
) const
{
    Info<< "Detecting near surfaces ..." << endl;
 
    const pointField& localPoints = pp.localPoints();
    const labelList& meshPoints = pp.meshPoints();
    const refinementSurfaces& surfaces = meshRefiner_.surfaces();
    const fvMesh& mesh = meshRefiner_.mesh();
 
    //const scalarField edgeLen(calcEdgeLen(pp));
    //
    //
    //{
    //    const scalar cos45 = Foam::cos(degToRad(45));
    //    vector n(cos45, cos45, cos45);
    //    n /= mag(n);
    //
    //    pointField start(14*pp.nPoints());
    //    pointField end(start.size());
    //
    //    label rayI = 0;
    //    forAll(localPoints, pointI)
    //    {
    //        const point& pt = localPoints[pointI];
    //
    //        // Along coordinate axes
    //
    //        {
    //            start[rayI] = pt;
    //            point& endPt = end[rayI++];
    //            endPt = pt;
    //            endPt.x() -= edgeLen[pointI];
    //        }
    //        {
    //            start[rayI] = pt;
    //            point& endPt = end[rayI++];
    //            endPt = pt;
    //            endPt.x() += edgeLen[pointI];
    //        }
    //        {
    //            start[rayI] = pt;
    //            point& endPt = end[rayI++];
    //            endPt = pt;
    //            endPt.y() -= edgeLen[pointI];
    //        }
    //        {
    //            start[rayI] = pt;
    //            point& endPt = end[rayI++];
    //            endPt = pt;
    //            endPt.y() += edgeLen[pointI];
    //        }
    //        {
    //            start[rayI] = pt;
    //            point& endPt = end[rayI++];
    //            endPt = pt;
    //            endPt.z() -= edgeLen[pointI];
    //        }
    //        {
    //            start[rayI] = pt;
    //            point& endPt = end[rayI++];
    //            endPt = pt;
    //            endPt.z() += edgeLen[pointI];
    //        }
    //
    //        // At 45 degrees
    //
    //        const vector vec(edgeLen[pointI]*n);
    //
    //        {
    //            start[rayI] = pt;
    //            point& endPt = end[rayI++];
    //            endPt = pt;
    //            endPt.x() += vec.x();
    //            endPt.y() += vec.y();
    //            endPt.z() += vec.z();
    //        }
    //        {
    //            start[rayI] = pt;
    //            point& endPt = end[rayI++];
    //            endPt = pt;
    //            endPt.x() -= vec.x();
    //            endPt.y() += vec.y();
    //            endPt.z() += vec.z();
    //        }
    //        {
    //            start[rayI] = pt;
    //            point& endPt = end[rayI++];
    //            endPt = pt;
    //            endPt.x() += vec.x();
    //            endPt.y() -= vec.y();
    //            endPt.z() += vec.z();
    //        }
    //        {
    //            start[rayI] = pt;
    //            point& endPt = end[rayI++];
    //            endPt = pt;
    //            endPt.x() -= vec.x();
    //            endPt.y() -= vec.y();
    //            endPt.z() += vec.z();
    //        }
    //        {
    //            start[rayI] = pt;
    //            point& endPt = end[rayI++];
    //            endPt = pt;
    //            endPt.x() += vec.x();
    //            endPt.y() += vec.y();
    //            endPt.z() -= vec.z();
    //        }
    //        {
    //            start[rayI] = pt;
    //            point& endPt = end[rayI++];
    //            endPt = pt;
    //            endPt.x() -= vec.x();
    //            endPt.y() += vec.y();
    //            endPt.z() -= vec.z();
    //        }
    //        {
    //            start[rayI] = pt;
    //            point& endPt = end[rayI++];
    //            endPt = pt;
    //            endPt.x() += vec.x();
    //            endPt.y() -= vec.y();
    //            endPt.z() -= vec.z();
    //        }
    //        {
    //            start[rayI] = pt;
    //            point& endPt = end[rayI++];
    //            endPt = pt;
    //            endPt.x() -= vec.x();
    //            endPt.y() -= vec.y();
    //            endPt.z() -= vec.z();
    //        }
    //    }
    //
    //    labelList surface1;
    //    List<pointIndexHit> hit1;
    //    labelList region1;
    //    vectorField normal1;
    //
    //    labelList surface2;
    //    List<pointIndexHit> hit2;
    //    labelList region2;
    //    vectorField normal2;
    //    surfaces.findNearestIntersection
    //    (
    //        unzonedSurfaces,    // surfacesToTest,
    //        start,
    //        end,
    //
    //        surface1,
    //        hit1,
    //        region1,
    //        normal1,
    //
    //        surface2,
    //        hit2,
    //        region2,
    //        normal2
    //    );
    //
    //    // All intersections
    //    {
    //        OBJstream str
    //        (
    //            mesh.time().path()
    //          / "surfaceHits_" + meshRefiner_.timeName() + ".obj"
    //        );
    //
    //        Info<< "Dumping intersections with rays to " << str.name()
    //            << endl;
    //
    //        forAll(hit1, i)
    //        {
    //            if (hit1[i].hit())
    //            {
    //                str.write(linePointRef(start[i], hit1[i].hitPoint()));
    //            }
    //            if (hit2[i].hit())
    //            {
    //                str.write(linePointRef(start[i], hit2[i].hitPoint()));
    //            }
    //        }
    //    }
    //
    //    // Co-planar intersections
    //    {
    //        OBJstream str
    //        (
    //            mesh.time().path()
    //          / "coplanarHits_" + meshRefiner_.timeName() + ".obj"
    //        );
    //
    //        Info<< "Dumping intersections with co-planar surfaces to "
    //            << str.name() << endl;
    //
    //        forAll(localPoints, pointI)
    //        {
    //            bool hasNormal = false;
    //            point surfPointA;
    //            vector surfNormalA;
    //            point surfPointB;
    //            vector surfNormalB;
    //
    //            bool isCoplanar = false;
    //
    //            label rayI = 14*pointI;
    //            for (label i = 0; i < 14; i++)
    //            {
    //                if (hit1[rayI].hit())
    //                {
    //                    const point& pt = hit1[rayI].hitPoint();
    //                    const vector& n = normal1[rayI];
    //
    //                    if (!hasNormal)
    //                    {
    //                        hasNormal = true;
    //                        surfPointA = pt;
    //                        surfNormalA = n;
    //                    }
    //                    else
    //                    {
    //                        if
    //                        (
    //                            meshRefiner_.isGap
    //                            (
    //                                planarCos,
    //                                surfPointA,
    //                                surfNormalA,
    //                                pt,
    //                                n
    //                            )
    //                        )
    //                        {
    //                            isCoplanar = true;
    //                            surfPointB = pt;
    //                            surfNormalB = n;
    //                            break;
    //                        }
    //                    }
    //                }
    //                if (hit2[rayI].hit())
    //                {
    //                    const point& pt = hit2[rayI].hitPoint();
    //                    const vector& n = normal2[rayI];
    //
    //                    if (!hasNormal)
    //                    {
    //                        hasNormal = true;
    //                        surfPointA = pt;
    //                        surfNormalA = n;
    //                    }
    //                    else
    //                    {
    //                        if
    //                        (
    //                            meshRefiner_.isGap
    //                            (
    //                                planarCos,
    //                                surfPointA,
    //                                surfNormalA,
    //                                pt,
    //                                n
    //                            )
    //                        )
    //                        {
    //                            isCoplanar = true;
    //                            surfPointB = pt;
    //                            surfNormalB = n;
    //                            break;
    //                        }
    //                    }
    //                }
    //
    //                rayI++;
    //            }
    //
    //            if (isCoplanar)
    //            {
    //                str.write(linePointRef(surfPointA, surfPointB));
    //            }
    //        }
    //    }
    //}
 
 
    const pointField avgCc(avgCellCentres(mesh, pp));
 
    // Construct rays through localPoints to beyond cell centre
    pointField start(pp.nPoints());
    pointField end(pp.nPoints());
    forAll(localPoints, pointI)
    {
        const point& pt = localPoints[pointI];
        const vector d = 2*(avgCc[pointI]-pt);
        start[pointI] = pt - d;
        end[pointI] = pt + d;
    }
 
 
    autoPtr<OBJstream> gapStr;
    if (debug&meshRefinement::ATTRACTION)
    {
        gapStr.reset
        (
            new OBJstream
            (
                mesh.time().path()
              / "detectNearSurfaces_" + meshRefiner_.timeName() + ".obj"
            )
        );
    }
 
 
    const PackedBoolList isPatchMasterPoint
    (
        meshRefinement::getMasterPoints
        (
            mesh,
            meshPoints
        )
    );
 
    label nOverride = 0;
 
    // 1. All points to non-interface surfaces
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    {
        const labelList unzonedSurfaces =
            surfaceZonesInfo::getUnnamedSurfaces
            (
                meshRefiner_.surfaces().surfZones()
            );
 
        // Do intersection test
        labelList surface1;
        List<pointIndexHit> hit1;
        labelList region1;
        vectorField normal1;
 
        labelList surface2;
        List<pointIndexHit> hit2;
        labelList region2;
        vectorField normal2;
        surfaces.findNearestIntersection
        (
            unzonedSurfaces,
            start,
            end,
 
            surface1,
            hit1,
            region1,
            normal1,
 
            surface2,
            hit2,
            region2,
            normal2
        );
 
 
        forAll(localPoints, pointI)
        {
            // Current location
            const point& pt = localPoints[pointI];
 
            bool override = false;
 
            //if (hit1[pointI].hit())
            //{
            //    if
            //    (
            //        meshRefiner_.isGap
            //        (
            //            planarCos,
            //            nearestPoint[pointI],
            //            nearestNormal[pointI],
            //            hit1[pointI].hitPoint(),
            //            normal1[pointI]
            //        )
            //    )
            //    {
            //        disp[pointI] = hit1[pointI].hitPoint()-pt;
            //        override = true;
            //    }
            //}
            //if (hit2[pointI].hit())
            //{
            //    if
            //    (
            //        meshRefiner_.isGap
            //        (
            //            planarCos,
            //            nearestPoint[pointI],
            //            nearestNormal[pointI],
            //            hit2[pointI].hitPoint(),
            //            normal2[pointI]
            //        )
            //    )
            //    {
            //        disp[pointI] = hit2[pointI].hitPoint()-pt;
            //        override = true;
            //    }
            //}
 
            if (hit1[pointI].hit() && hit2[pointI].hit())
            {
                if
                (
                    meshRefiner_.isGap
                    (
                        planarCos,
                        hit1[pointI].hitPoint(),
                        normal1[pointI],
                        hit2[pointI].hitPoint(),
                        normal2[pointI]
                    )
                )
                {
                    // TBD: check if the attraction (to nearest) would attract
                    // good enough and not override attraction
 
                    if (gapStr.valid())
                    {
                        const point& intPt = hit2[pointI].hitPoint();
                        gapStr().write(linePointRef(pt, intPt));
                    }
 
                    // Choose hit2 : nearest to end point (so inside the domain)
                    disp[pointI] = hit2[pointI].hitPoint()-pt;
                    override = true;
                }
            }
 
            if (override && isPatchMasterPoint[pointI])
            {
                nOverride++;
            }
        }
    }
 
 
    // 2. All points on zones to their respective surface
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
    {
        // Surfaces with zone information
        const PtrList<surfaceZonesInfo>& surfZones = surfaces.surfZones();
 
        const labelList zonedSurfaces = surfaceZonesInfo::getNamedSurfaces
        (
            surfZones
        );
 
        forAll(zonedSurfaces, i)
        {
            label zoneSurfI = zonedSurfaces[i];
 
            const word& faceZoneName = surfZones[zoneSurfI].faceZoneName();
 
            const labelList surfacesToTest(1, zoneSurfI);
 
            // Get indices of points both on faceZone and on pp.
            labelList zonePointIndices
            (
                getZoneSurfacePoints
                (
                    mesh,
                    pp,
                    faceZoneName
                )
            );
 
            // Do intersection test
            labelList surface1;
            List<pointIndexHit> hit1;
            labelList region1;
            vectorField normal1;
 
            labelList surface2;
            List<pointIndexHit> hit2;
            labelList region2;
            vectorField normal2;
            surfaces.findNearestIntersection
            (
                surfacesToTest,
                pointField(start, zonePointIndices),
                pointField(end, zonePointIndices),
 
                surface1,
                hit1,
                region1,
                normal1,
 
                surface2,
                hit2,
                region2,
                normal2
            );
 
 
            forAll(hit1, i)
            {
                label pointI = zonePointIndices[i];
 
                // Current location
                const point& pt = localPoints[pointI];
 
                bool override = false;
 
                //if (hit1[i].hit())
                //{
                //    if
                //    (
                //        meshRefiner_.isGap
                //        (
                //            planarCos,
                //            nearestPoint[pointI],
                //            nearestNormal[pointI],
                //            hit1[i].hitPoint(),
                //            normal1[i]
                //        )
                //    )
                //    {
                //        disp[pointI] = hit1[i].hitPoint()-pt;
                //        override = true;
                //    }
                //}
                //if (hit2[i].hit())
                //{
                //    if
                //    (
                //        meshRefiner_.isGap
                //        (
                //            planarCos,
                //            nearestPoint[pointI],
                //            nearestNormal[pointI],
                //            hit2[i].hitPoint(),
                //            normal2[i]
                //        )
                //    )
                //    {
                //        disp[pointI] = hit2[i].hitPoint()-pt;
                //        override = true;
                //    }
                //}
 
                if (hit1[i].hit() && hit2[i].hit())
                {
                    if
                    (
                        meshRefiner_.isGap
                        (
                            planarCos,
                            hit1[i].hitPoint(),
                            normal1[i],
                            hit2[i].hitPoint(),
                            normal2[i]
                        )
                    )
                    {
                        if (gapStr.valid())
                        {
                            const point& intPt = hit2[i].hitPoint();
                            gapStr().write(linePointRef(pt, intPt));
                        }
 
                        disp[pointI] = hit2[i].hitPoint()-pt;
                        override = true;
                    }
                }
 
                if (override && isPatchMasterPoint[pointI])
                {
                    nOverride++;
                }
            }
        }
    }
 
    Info<< "Overriding nearest with intersection of close gaps at "
        << returnReduce(nOverride, sumOp<label>())
        << " out of " << returnReduce(pp.nPoints(), sumOp<label>())
        << " points." << endl;
}
 
 
void Foam::autoSnapDriver::calcNearestSurface
(
    const refinementSurfaces& surfaces,
 
    const labelList& surfacesToTest,
    const labelListList& regionsToTest,
 
    const pointField& localPoints,
    const labelList& zonePointIndices,
 
    scalarField& minSnapDist,
    labelList& snapSurf,
    vectorField& patchDisp,
 
    // Optional: nearest point, normal
    pointField& nearestPoint,
    vectorField& nearestNormal
)
{
    // Find nearest for points both on faceZone and pp.
    List<pointIndexHit> hitInfo;
    labelList hitSurface;
 
    if (nearestNormal.size() == localPoints.size())
    {
        labelList hitRegion;
        vectorField hitNormal;
        surfaces.findNearestRegion
        (
            surfacesToTest,
            regionsToTest,
 
            pointField(localPoints, zonePointIndices),
            sqr(scalarField(minSnapDist, zonePointIndices)),
 
            hitSurface,
            hitInfo,
            hitRegion,
            hitNormal
        );
 
        forAll(hitInfo, i)
        {
            if (hitInfo[i].hit())
            {
                label pointI = zonePointIndices[i];
                nearestPoint[pointI] = hitInfo[i].hitPoint();
                nearestNormal[pointI] = hitNormal[i];
            }
        }
    }
    else
    {
        surfaces.findNearest
        (
            surfacesToTest,
            regionsToTest,
 
            pointField(localPoints, zonePointIndices),
            sqr(scalarField(minSnapDist, zonePointIndices)),
 
            hitSurface,
            hitInfo
        );
    }
 
    forAll(hitInfo, i)
    {
        if (hitInfo[i].hit())
        {
            label pointI = zonePointIndices[i];
 
            patchDisp[pointI] = hitInfo[i].hitPoint() - localPoints[pointI];
            minSnapDist[pointI] = mag(patchDisp[pointI]);
            snapSurf[pointI] = hitSurface[i];
        }
    }
}
 
 
Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
(
    const bool strictRegionSnap,
    const meshRefinement& meshRefiner,
    const labelList& globalToMasterPatch,
    const labelList& globalToSlavePatch,
    const scalarField& snapDist,
    const indirectPrimitivePatch& pp,
    pointField& nearestPoint,
    vectorField& nearestNormal
)
{
    Info<< "Calculating patchDisplacement as distance to nearest surface"
        << " point ..." << endl;
    if (strictRegionSnap)
    {
        Info<< "    non-zone points : attract to local region on surface only"
            << nl
            << "    zone points     : attract to local region on surface only"
            << nl
            << endl;
    }
    else
    {
        Info<< "    non-zone points :"
            << " attract to nearest of all non-zone surfaces"
            << nl
            << "    zone points     : attract to zone surface only" << nl
            << endl;
    }
 
 
    const pointField& localPoints = pp.localPoints();
    const refinementSurfaces& surfaces = meshRefiner.surfaces();
    const fvMesh& mesh = meshRefiner.mesh();
 
    // Displacement per patch point
    vectorField patchDisp(localPoints.size(), vector::zero);
 
    if (returnReduce(localPoints.size(), sumOp<label>()) > 0)
    {
        // Current surface snapped to. Used to check whether points have been
        // snapped at all
        labelList snapSurf(localPoints.size(), -1);
 
        // Current best snap distance (since point might be on multiple
        // regions)
        scalarField minSnapDist(snapDist);
 
 
        if (strictRegionSnap)
        {
            // Attract patch points to same region only
 
            forAll(surfaces.surfaces(), surfI)
            {
                label geomI = surfaces.surfaces()[surfI];
                label nRegions = surfaces.geometry()[geomI].regions().size();
 
                const labelList surfacesToTest(1, surfI);
 
                for (label regionI = 0; regionI < nRegions; regionI++)
                {
                    label globalI = surfaces.globalRegion(surfI, regionI);
                    label masterPatchI = globalToMasterPatch[globalI];
 
                    // Get indices of points both on patch and on pp
                    labelList zonePointIndices
                    (
                        getFacePoints
                        (
                            pp,
                            mesh.boundaryMesh()[masterPatchI]
                        )
                    );
 
                    calcNearestSurface
                    (
                        surfaces,
 
                        surfacesToTest,
                        labelListList(1, labelList(1, regionI)), //regionsToTest
 
                        localPoints,
                        zonePointIndices,
 
                        minSnapDist,
                        snapSurf,
                        patchDisp,
 
                        // Optional: nearest point, normal
                        nearestPoint,
                        nearestNormal
                    );
 
                    if (globalToSlavePatch[globalI] != masterPatchI)
                    {
                        label slavePatchI = globalToSlavePatch[globalI];
 
                        // Get indices of points both on patch and on pp
                        labelList zonePointIndices
                        (
                            getFacePoints
                            (
                                pp,
                                mesh.boundaryMesh()[slavePatchI]
                            )
                        );
 
                        calcNearestSurface
                        (
                            surfaces,
 
                            surfacesToTest,
                            labelListList(1, labelList(1, regionI)),
 
                            localPoints,
                            zonePointIndices,
 
                            minSnapDist,
                            snapSurf,
                            patchDisp,
 
                            // Optional: nearest point, normal
                            nearestPoint,
                            nearestNormal
                        );
                    }
                }
            }
        }
        else
        {
            // Divide surfaces into zoned and unzoned
            const labelList unzonedSurfaces =
                surfaceZonesInfo::getUnnamedSurfaces
                (
                    meshRefiner.surfaces().surfZones()
                );
 
 
            // 1. All points to non-interface surfaces
            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
            List<pointIndexHit> hitInfo;
            labelList hitSurface;
 
            if (nearestNormal.size() == localPoints.size())
            {
                labelList hitRegion;
                vectorField hitNormal;
                surfaces.findNearestRegion
                (
                    unzonedSurfaces,
                    localPoints,
                    sqr(snapDist),
                    hitSurface,
                    hitInfo,
                    hitRegion,
                    hitNormal
                );
 
                forAll(hitInfo, pointI)
                {
                    if (hitInfo[pointI].hit())
                    {
                        nearestPoint[pointI] = hitInfo[pointI].hitPoint();
                        nearestNormal[pointI] = hitNormal[pointI];
                    }
                }
            }
            else
            {
                surfaces.findNearest
                (
                    unzonedSurfaces,
                    localPoints,
                    sqr(snapDist),        // sqr of attract distance
                    hitSurface,
                    hitInfo
                );
            }
 
            forAll(hitInfo, pointI)
            {
                if (hitInfo[pointI].hit())
                {
                    patchDisp[pointI] =
                        hitInfo[pointI].hitPoint()
                      - localPoints[pointI];
 
                    snapSurf[pointI] = hitSurface[pointI];
                }
            }
 
 
            const labelList zonedSurfaces =
                surfaceZonesInfo::getNamedSurfaces
                (
                    meshRefiner.surfaces().surfZones()
                );
 
 
            // 2. All points on zones to their respective surface
            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
            // Surfaces with zone information
            const PtrList<surfaceZonesInfo>& surfZones = surfaces.surfZones();
 
            forAll(zonedSurfaces, i)
            {
                label surfI = zonedSurfaces[i];
 
                const word& faceZoneName = surfZones[surfI].faceZoneName();
 
                const labelList surfacesToTest(1, surfI);
 
                label geomI = surfaces.surfaces()[surfI];
                label nRegions = surfaces.geometry()[geomI].regions().size();
 
 
                // Get indices of points both on faceZone and on pp.
                labelList zonePointIndices
                (
                    getZoneSurfacePoints
                    (
                        mesh,
                        pp,
                        faceZoneName
                    )
                );
 
 
                calcNearestSurface
                (
                    surfaces,
 
                    surfacesToTest,
                    labelListList(1, identity(nRegions)),
 
                    localPoints,
                    zonePointIndices,
 
                    minSnapDist,
                    snapSurf,
                    patchDisp,
 
                    // Optional: nearest point, normal
                    nearestPoint,
                    nearestNormal
                );
            }
        }
 
 
        // Check if all points are being snapped
        forAll(snapSurf, pointI)
        {
            if (snapSurf[pointI] == -1)
            {
                WarningInFunction
                    << "For point:" << pointI
                    << " coordinate:" << localPoints[pointI]
                    << " did not find any surface within:"
                    << minSnapDist[pointI]
                    << " metre." << endl;
            }
        }
 
        {
            const PackedBoolList isPatchMasterPoint
            (
                meshRefinement::getMasterPoints
                (
                    mesh,
                    pp.meshPoints()
                )
            );
 
            scalarField magDisp(mag(patchDisp));
 
            Info<< "Wanted displacement : average:"
                <<  meshRefinement::gAverage(isPatchMasterPoint, magDisp)
                << " min:" << gMin(magDisp)
                << " max:" << gMax(magDisp) << endl;
        }
    }
 
    Info<< "Calculated surface displacement in = "
        << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
 
 
    // Limit amount of movement. Can not happen for triSurfaceMesh but
    // can happen for some analytical shapes?
    forAll(patchDisp, patchPointI)
    {
        scalar magDisp = mag(patchDisp[patchPointI]);
 
        if (magDisp > snapDist[patchPointI])
        {
            patchDisp[patchPointI] *= snapDist[patchPointI] / magDisp;
 
            Pout<< "Limiting displacement for " << patchPointI
                << " from " << magDisp << " to " << snapDist[patchPointI]
                << endl;
        }
    }
 
    // Points on zones in one domain but only present as point on other
    // will not do condition 2 on all. Sync explicitly.
    syncTools::syncPointList
    (
        mesh,
        pp.meshPoints(),
        patchDisp,
        minMagSqrEqOp<point>(),         // combine op
        vector(GREAT, GREAT, GREAT)     // null value (note: cannot use VGREAT)
    );
 
    return patchDisp;
}
 
 
void Foam::autoSnapDriver::smoothDisplacement
(
    const snapParameters& snapParams,
    motionSmoother& meshMover
) const
{
    const fvMesh& mesh = meshRefiner_.mesh();
    const indirectPrimitivePatch& pp = meshMover.patch();
 
    Info<< "Smoothing displacement ..." << endl;
 
    // Set edge diffusivity as inverse of distance to patch
    scalarField edgeGamma(1.0/(edgePatchDist(meshMover.pMesh(), pp) + SMALL));
    //scalarField edgeGamma(mesh.nEdges(), 1.0);
    //scalarField edgeGamma(wallGamma(mesh, pp, 10, 1));
 
    // Get displacement field
    pointVectorField& disp = meshMover.displacement();
 
    for (label iter = 0; iter < snapParams.nSmoothDispl(); iter++)
    {
        if ((iter % 10) == 0)
        {
            Info<< "Iteration " << iter << endl;
        }
        pointVectorField oldDisp(disp);
        meshMover.smooth(oldDisp, edgeGamma, disp);
    }
    Info<< "Displacement smoothed in = "
        << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
 
    if (debug&meshRefinement::MESH)
    {
        const_cast<Time&>(mesh.time())++;
        Info<< "Writing smoothed mesh to time " << meshRefiner_.timeName()
            << endl;
 
        // Moving mesh creates meshPhi. Can be cleared out by a mesh.clearOut
        // but this will also delete all pointMesh but not pointFields which
        // gives an illegal situation.
 
        meshRefiner_.write
        (
            meshRefinement::debugType(debug),
            meshRefinement::writeType
            (
                meshRefinement::writeLevel()
              | meshRefinement::WRITEMESH
            ),
            mesh.time().path()/meshRefiner_.timeName()
        );
        Info<< "Writing displacement field ..." << endl;
        disp.write();
        tmp<pointScalarField> magDisp(mag(disp));
        magDisp().write();
 
        Info<< "Writing actual patch displacement ..." << endl;
        vectorField actualPatchDisp(disp, pp.meshPoints());
        dumpMove
        (
            mesh.time().path()
          / "actualPatchDisplacement_" + meshRefiner_.timeName() + ".obj",
            pp.localPoints(),
            pp.localPoints() + actualPatchDisp
        );
    }
}
 
 
bool Foam::autoSnapDriver::scaleMesh
(
    const snapParameters& snapParams,
    const label nInitErrors,
    const List<labelPair>& baffles,
    motionSmoother& meshMover
)
{
    const fvMesh& mesh = meshRefiner_.mesh();
 
    // Relax displacement until correct mesh
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    labelList checkFaces(identity(mesh.nFaces()));
 
    scalar oldErrorReduction = -1;
 
    bool meshOk = false;
 
    Info<< "Moving mesh ..." << endl;
    for (label iter = 0; iter < 2*snapParams.nSnap(); iter++)
    {
        Info<< nl << "Iteration " << iter << endl;
 
        if (iter == snapParams.nSnap())
        {
            Info<< "Displacement scaling for error reduction set to 0." << endl;
            oldErrorReduction = meshMover.setErrorReduction(0.0);
        }
 
        meshOk = meshMover.scaleMesh(checkFaces, baffles, true, nInitErrors);
 
        if (meshOk)
        {
            Info<< "Successfully moved mesh" << endl;
            break;
        }
        if (debug&meshRefinement::MESH)
        {
            const_cast<Time&>(mesh.time())++;
            Info<< "Writing scaled mesh to time " << meshRefiner_.timeName()
                << endl;
            mesh.write();
 
            Info<< "Writing displacement field ..." << endl;
            meshMover.displacement().write();
            tmp<pointScalarField> magDisp(mag(meshMover.displacement()));
            magDisp().write();
        }
    }
 
    if (oldErrorReduction >= 0)
    {
        meshMover.setErrorReduction(oldErrorReduction);
    }
    Info<< "Moved mesh in = "
        << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
 
    return meshOk;
}
 
 
// After snapping: correct patching according to nearest surface.
// Code is very similar to calcNearestSurface.
// - calculate face-wise snap distance as max of point-wise
// - calculate face-wise nearest surface point
// - repatch face according to patch for surface point.
Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::repatchToSurface
(
    const snapParameters& snapParams,
    const labelList& adaptPatchIDs,
    const labelList& preserveFaces
)
{
    const fvMesh& mesh = meshRefiner_.mesh();
    const refinementSurfaces& surfaces = meshRefiner_.surfaces();
 
    Info<< "Repatching faces according to nearest surface ..." << endl;
 
    // Get the labels of added patches.
    autoPtr<indirectPrimitivePatch> ppPtr
    (
        meshRefinement::makePatch
        (
            mesh,
            adaptPatchIDs
        )
    );
    indirectPrimitivePatch& pp = ppPtr();
 
    // Divide surfaces into zoned and unzoned
    labelList zonedSurfaces =
        surfaceZonesInfo::getNamedSurfaces(surfaces.surfZones());
    labelList unzonedSurfaces =
        surfaceZonesInfo::getUnnamedSurfaces(surfaces.surfZones());
 
 
    // Faces that do not move
    PackedBoolList isZonedFace(mesh.nFaces());
    {
        // 1. Preserve faces in preserveFaces list
        forAll(preserveFaces, faceI)
        {
            if (preserveFaces[faceI] != -1)
            {
                isZonedFace.set(faceI, 1);
            }
        }
 
        // 2. All faces on zoned surfaces
        const PtrList<surfaceZonesInfo>& surfZones = surfaces.surfZones();
        const faceZoneMesh& fZones = mesh.faceZones();
 
        forAll(zonedSurfaces, i)
        {
            const label zoneSurfI = zonedSurfaces[i];
            const faceZone& fZone = fZones[surfZones[zoneSurfI].faceZoneName()];
 
            forAll(fZone, i)
            {
                isZonedFace.set(fZone[i], 1);
            }
        }
    }
 
 
    // Determine per pp face which patch it should be in
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
    // Patch that face should be in
    labelList closestPatch(pp.size(), -1);
    {
        // face snap distance as max of point snap distance
        scalarField faceSnapDist(pp.size(), -GREAT);
        {
            // Distance to attract to nearest feature on surface
            const scalarField snapDist
            (
                calcSnapDistance
                (
                    mesh,
                    snapParams,
                    pp
                )
            );
 
            const faceList& localFaces = pp.localFaces();
 
            forAll(localFaces, faceI)
            {
                const face& f = localFaces[faceI];
 
                forAll(f, fp)
                {
                    faceSnapDist[faceI] = max
                    (
                        faceSnapDist[faceI],
                        snapDist[f[fp]]
                    );
                }
            }
        }
 
        pointField localFaceCentres(mesh.faceCentres(), pp.addressing());
 
        // Get nearest surface and region
        labelList hitSurface;
        labelList hitRegion;
        surfaces.findNearestRegion
        (
            unzonedSurfaces,
            localFaceCentres,
            sqr(faceSnapDist),    // sqr of attract distance
            hitSurface,
            hitRegion
        );
 
        // Get patch
        forAll(pp, i)
        {
            label faceI = pp.addressing()[i];
 
            if (hitSurface[i] != -1 && !isZonedFace.get(faceI))
            {
                closestPatch[i] = globalToMasterPatch_
                [
                    surfaces.globalRegion
                    (
                        hitSurface[i],
                        hitRegion[i]
                    )
                ];
            }
        }
    }
 
 
    // Change those faces for which there is a different closest patch
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
    labelList ownPatch(mesh.nFaces(), -1);
    labelList neiPatch(mesh.nFaces(), -1);
 
    const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
    forAll(patches, patchI)
    {
        const polyPatch& pp = patches[patchI];
 
        forAll(pp, i)
        {
            ownPatch[pp.start()+i] = patchI;
            neiPatch[pp.start()+i] = patchI;
        }
    }
 
    label nChanged = 0;
    forAll(closestPatch, i)
    {
        label faceI = pp.addressing()[i];
 
        if (closestPatch[i] != -1 && closestPatch[i] != ownPatch[faceI])
        {
            ownPatch[faceI] = closestPatch[i];
            neiPatch[faceI] = closestPatch[i];
            nChanged++;
        }
    }
 
    Info<< "Repatched " << returnReduce(nChanged, sumOp<label>())
        << " faces in = " << mesh.time().cpuTimeIncrement() << " s\n" << nl
        << endl;
 
    return meshRefiner_.createBaffles(ownPatch, neiPatch);
}
 
 
void Foam::autoSnapDriver::detectWarpedFaces
(
    const scalar featureCos,
    const indirectPrimitivePatch& pp,
 
    DynamicList<label>& splitFaces,
    DynamicList<labelPair>& splits
) const
{
    const fvMesh& mesh = meshRefiner_.mesh();
    const faceList& localFaces = pp.localFaces();
    const pointField& localPoints = pp.localPoints();
    const labelList& bFaces = pp.addressing();
 
    splitFaces.clear();
    splitFaces.setCapacity(bFaces.size());
    splits.clear();
    splits.setCapacity(bFaces.size());
 
    // Determine parallel consistent normals on points
    const vectorField pointNormals(PatchTools::pointNormals(mesh, pp));
 
    face f0(4);
    face f1(4);
 
    forAll(localFaces, faceI)
    {
        const face& f = localFaces[faceI];
 
        if (f.size() >= 4)
        {
            // See if splitting face across diagonal would make two faces with
            // biggish normal angle
 
            labelPair minDiag(-1, -1);
            scalar minCos(GREAT);
 
            for (label startFp = 0; startFp < f.size()-2; startFp++)
            {
                label minFp = f.rcIndex(startFp);
 
                for
                (
                    label endFp = f.fcIndex(f.fcIndex(startFp));
                    endFp < f.size() && endFp != minFp;
                    endFp++
                )
                {
                    // Form two faces
                    f0.setSize(endFp-startFp+1);
                    label i0 = 0;
                    for (label fp = startFp; fp <= endFp; fp++)
                    {
                        f0[i0++] = f[fp];
                    }
                    f1.setSize(f.size()+2-f0.size());
                    label i1 = 0;
                    for (label fp = endFp; fp != startFp; fp = f.fcIndex(fp))
                    {
                        f1[i1++] = f[fp];
                    }
                    f1[i1++] = f[startFp];
 
                    //Info<< "Splitting face:" << f << " into f0:" << f0
                    //    << " f1:" << f1 << endl;
 
                    vector n0 = f0.normal(localPoints);
                    scalar n0Mag = mag(n0);
                    vector n1 = f1.normal(localPoints);
                    scalar n1Mag = mag(n1);
 
                    if (n0Mag > ROOTVSMALL && n1Mag > ROOTVSMALL)
                    {
                        scalar cosAngle = (n0/n0Mag) & (n1/n1Mag);
                        if (cosAngle < minCos)
                        {
                            minCos = cosAngle;
                            minDiag = labelPair(startFp, endFp);
                        }
                    }
                }
            }
 
 
            if (minCos < featureCos)
            {
                splitFaces.append(bFaces[faceI]);
                splits.append(minDiag);
            }
        }
    }
}
 
 
Foam::labelList Foam::autoSnapDriver::getInternalOrBaffleDuplicateFace() const
{
    const fvMesh& mesh = meshRefiner_.mesh();
 
    labelList internalOrBaffleFaceZones;
    {
        List<surfaceZonesInfo::faceZoneType> fzTypes(2);
        fzTypes[0] = surfaceZonesInfo::INTERNAL;
        fzTypes[1] = surfaceZonesInfo::BAFFLE;
        internalOrBaffleFaceZones = meshRefiner_.getZones(fzTypes);
    }
 
    List<labelPair> baffles
    (
        meshRefiner_.subsetBaffles
        (
            mesh,
            internalOrBaffleFaceZones,
            localPointRegion::findDuplicateFacePairs(mesh)
        )
    );
 
    labelList faceToDuplicate(mesh.nFaces(), -1);
    forAll(baffles, i)
    {
        const labelPair& p = baffles[i];
        faceToDuplicate[p[0]] = p[1];
        faceToDuplicate[p[1]] = p[0];
    }
 
    return faceToDuplicate;
}
 
 
void Foam::autoSnapDriver::doSnap
(
    const dictionary& snapDict,
    const dictionary& motionDict,
    const bool mergePatchFaces,
    const scalar featureCos,
    const scalar planarAngle,
    const snapParameters& snapParams
)
{
    fvMesh& mesh = meshRefiner_.mesh();
 
    Info<< nl
        << "Morphing phase" << nl
        << "--------------" << nl
        << endl;
 
    // faceZone handling
    // ~~~~~~~~~~~~~~~~~
    //
    // We convert all faceZones into baffles during snapping so we can use
    // a standard mesh motion (except for the mesh checking which for baffles
    // created from internal faces should check across the baffles). The state
    // is stored in two variables:
    //      baffles : pairs of boundary faces
    //      duplicateFace : from mesh face to its baffle colleague (or -1 for
    //                      normal faces)
    // There are three types of faceZones according to the faceType property:
    //
    // internal
    // --------
    // - baffles: need to be checked across
    // - duplicateFace: from face to duplicate face. Contains
    //   all faces on faceZone to prevents merging patch faces.
    //
    // baffle
    // ------
    // - baffles: no need to be checked across
    // - duplicateFace: contains all faces on faceZone to prevent
    //   merging patch faces.
    //
    // boundary
    // --------
    // - baffles: no need to be checked across. Also points get duplicated
    //            so will no longer be baffles
    // - duplicateFace: contains no faces on faceZone since both sides can
    //   merge faces independently.
 
 
 
    // faceZones of type internal
    const labelList internalFaceZones
    (
        meshRefiner_.getZones
        (
            List<surfaceZonesInfo::faceZoneType>
            (
                1,
                surfaceZonesInfo::INTERNAL
            )
        )
    );
 
 
    // Create baffles (pairs of faces that share the same points)
    // Baffles stored as owner and neighbour face that have been created.
    {
        List<labelPair> baffles;
        labelList originatingFaceZone;
        meshRefiner_.createZoneBaffles
        (
            identity(mesh.faceZones().size()),
            baffles,
            originatingFaceZone
        );
    }
 
    // Duplicate points on faceZones of type boundary
    meshRefiner_.dupNonManifoldBoundaryPoints();
 
 
    bool doFeatures = false;
    label nFeatIter = 1;
    if (snapParams.nFeatureSnap() > 0)
    {
        doFeatures = true;
        nFeatIter = snapParams.nFeatureSnap();
 
        Info<< "Snapping to features in " << nFeatIter
            << " iterations ..." << endl;
    }
 
 
    bool meshOk = false;
 
 
    // Get the labels of added patches.
    labelList adaptPatchIDs(meshRefiner_.meshedPatches());
 
 
 
    {
        autoPtr<indirectPrimitivePatch> ppPtr
        (
            meshRefinement::makePatch
            (
                mesh,
                adaptPatchIDs
            )
        );
 
 
        // Distance to attract to nearest feature on surface
        scalarField snapDist(calcSnapDistance(mesh, snapParams, ppPtr()));
 
 
        // Construct iterative mesh mover.
        Info<< "Constructing mesh displacer ..." << endl;
        Info<< "Using mesh parameters " << motionDict << nl << endl;
 
        autoPtr<motionSmoother> meshMoverPtr
        (
            new motionSmoother
            (
                mesh,
                ppPtr(),
                adaptPatchIDs,
                meshRefinement::makeDisplacementField
                (
                    pointMesh::New(mesh),
                    adaptPatchIDs
                ),
                motionDict
            )
        );
 
 
        // Check initial mesh
        Info<< "Checking initial mesh ..." << endl;
        labelHashSet wrongFaces(mesh.nFaces()/100);
        motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces);
        const label nInitErrors = returnReduce
        (
            wrongFaces.size(),
            sumOp<label>()
        );
 
        Info<< "Detected " << nInitErrors << " illegal faces"
            << " (concave, zero area or negative cell pyramid volume)"
            << endl;
 
 
        Info<< "Checked initial mesh in = "
            << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
 
        // Extract baffles across internal faceZones (for checking mesh quality
        // across
        labelPairList internalBaffles
        (
            meshRefiner_.subsetBaffles
            (
                mesh,
                internalFaceZones,
                localPointRegion::findDuplicateFacePairs(mesh)
            )
        );
 
 
 
        // Pre-smooth patch vertices (so before determining nearest)
        preSmoothPatch
        (
            meshRefiner_,
            snapParams,
            nInitErrors,
            internalBaffles,
            meshMoverPtr()
        );
 
 
 
        //- Only if in feature attraction mode:
        // Nearest feature
        vectorField patchAttraction;
        // Constraints at feature
        List<pointConstraint> patchConstraints;
 
 
        //- Any faces to split
        DynamicList<label> splitFaces;
        //- Indices in face to split across
        DynamicList<labelPair> splits;
 
 
        for (label iter = 0; iter < nFeatIter; iter++)
        {
            Info<< nl
                << "Morph iteration " << iter << nl
                << "-----------------" << endl;
 
            // Splitting iteration?
            bool doSplit = false;
            if
            (
                doFeatures
             && snapParams.nFaceSplitInterval() > 0
             && (
                    (iter == nFeatIter-1)
                 || (iter > 0 && (iter%snapParams.nFaceSplitInterval()) == 0)
                )
            )
            {
                doSplit = true;
            }
 
 
 
            indirectPrimitivePatch& pp = ppPtr();
            motionSmoother& meshMover = meshMoverPtr();
 
 
            // Calculate displacement at every patch point. Insert into
            // meshMover.
            // Calculate displacement at every patch point
            pointField nearestPoint;
            vectorField nearestNormal;
 
            if (snapParams.detectNearSurfacesSnap())
            {
                nearestPoint.setSize(pp.nPoints(), vector::max);
                nearestNormal.setSize(pp.nPoints(), vector::zero);
            }
 
            vectorField disp = calcNearestSurface
            (
                snapParams.strictRegionSnap(),  // attract points to region only
                meshRefiner_,
                globalToMasterPatch_,           // for if strictRegionSnap
                globalToSlavePatch_,            // for if strictRegionSnap
                snapDist,
                pp,
 
                nearestPoint,
                nearestNormal
            );
 
 
            // Override displacement at thin gaps
            if (snapParams.detectNearSurfacesSnap())
            {
                detectNearSurfaces
                (
                    Foam::cos(degToRad(planarAngle)),// planar cos for gaps
                    pp,
                    nearestPoint,   // surfacepoint from nearest test
                    nearestNormal,  // surfacenormal from nearest test
 
                    disp
                );
            }
 
            // Override displacement with feature edge attempt
            if (doFeatures)
            {
                splitFaces.clear();
                splits.clear();
                disp = calcNearestSurfaceFeature
                (
                    snapParams,
                    !doSplit,       // alignMeshEdges
                    iter,
                    featureCos,
                    scalar(iter+1)/nFeatIter,
 
                    snapDist,
                    disp,
                    nearestNormal,
                    meshMover,
 
                    patchAttraction,
                    patchConstraints,
 
                    splitFaces,
                    splits
                );
            }
 
            // Check for displacement being outwards.
            outwardsDisplacement(pp, disp);
 
            // Set initial distribution of displacement field (on patches)
            // from patchDisp and make displacement consistent with b.c.
            // on displacement pointVectorField.
            meshMover.setDisplacement(disp);
 
 
            if (debug&meshRefinement::ATTRACTION)
            {
                dumpMove
                (
                    mesh.time().path()
                  / "patchDisplacement_" + name(iter) + ".obj",
                    pp.localPoints(),
                    pp.localPoints() + disp
                );
            }
 
            // Get smoothly varying internal displacement field.
            smoothDisplacement(snapParams, meshMover);
 
            // Apply internal displacement to mesh.
            meshOk = scaleMesh
            (
                snapParams,
                nInitErrors,
                internalBaffles,
                meshMover
            );
 
            if (!meshOk)
            {
                WarningInFunction
                    << "Did not succesfully snap mesh."
                    << " Continuing to snap to resolve easy" << nl
                    << "    surfaces but the"
                    << " resulting mesh will not satisfy your quality"
                    << " constraints" << nl << endl;
            }
 
            if (debug&meshRefinement::MESH)
            {
                const_cast<Time&>(mesh.time())++;
                Info<< "Writing scaled mesh to time "
                    << meshRefiner_.timeName() << endl;
                meshRefiner_.write
                (
                    meshRefinement::debugType(debug),
                    meshRefinement::writeType
                    (
                        meshRefinement::writeLevel()
                      | meshRefinement::WRITEMESH
                    ),
                    mesh.time().path()/meshRefiner_.timeName()
                );
                Info<< "Writing displacement field ..." << endl;
                meshMover.displacement().write();
                tmp<pointScalarField> magDisp
                (
                    mag(meshMover.displacement())
                );
                magDisp().write();
            }
 
            // Use current mesh as base mesh
            meshMover.correct();
 
 
 
            // See if any faces need splitting
            label nTotalSplit = returnReduce(splitFaces.size(), sumOp<label>());
            if (nTotalSplit && doSplit)
            {
                // Filter out baffle faces from faceZones of type
                // internal/baffle
 
                labelList duplicateFace(getInternalOrBaffleDuplicateFace());
 
                {
                    labelList oldSplitFaces(splitFaces.xfer());
                    List<labelPair> oldSplits(splits.xfer());
                    forAll(oldSplitFaces, i)
                    {
                        if (duplicateFace[oldSplitFaces[i]] == -1)
                        {
                            splitFaces.append(oldSplitFaces[i]);
                            splits.append(oldSplits[i]);
                        }
                    }
                    nTotalSplit = returnReduce
                    (
                        splitFaces.size(),
                        sumOp<label>()
                    );
                }
 
                // Update mesh
                meshRefiner_.splitFacesUndo
                (
                    splitFaces,
                    splits,
                    motionDict,
 
                    duplicateFace,
                    internalBaffles
                );
 
                // Redo meshMover
                meshMoverPtr.clear();
                ppPtr.clear();
 
                // Update mesh mover
                ppPtr = meshRefinement::makePatch(mesh, adaptPatchIDs);
                meshMoverPtr.reset
                (
                    new motionSmoother
                    (
                        mesh,
                        ppPtr(),
                        adaptPatchIDs,
                        meshRefinement::makeDisplacementField
                        (
                            pointMesh::New(mesh),
                            adaptPatchIDs
                        ),
                        motionDict
                    )
                );
 
                // Update snapping distance
                snapDist = calcSnapDistance(mesh, snapParams, ppPtr());
 
 
                if (debug&meshRefinement::MESH)
                {
                    const_cast<Time&>(mesh.time())++;
                    Info<< "Writing split-faces mesh to time "
                        << meshRefiner_.timeName() << endl;
                    meshRefiner_.write
                    (
                        meshRefinement::debugType(debug),
                        meshRefinement::writeType
                        (
                            meshRefinement::writeLevel()
                          | meshRefinement::WRITEMESH
                        ),
                        mesh.time().path()/meshRefiner_.timeName()
                    );
                }
            }
 
 
            if (debug&meshRefinement::MESH)
            {
                forAll(internalBaffles, i)
                {
                    const labelPair& p = internalBaffles[i];
                    const point& fc0 = mesh.faceCentres()[p[0]];
                    const point& fc1 = mesh.faceCentres()[p[1]];
 
                    if (mag(fc0-fc1) > meshRefiner_.mergeDistance())
                    {
                        FatalErrorInFunction
                            << "Separated baffles : f0:" << p[0]
                            << " centre:" << fc0
                            << " f1:" << p[1] << " centre:" << fc1
                            << " distance:" << mag(fc0-fc1)
                            << exit(FatalError);
                    }
                }
            }
        }
    }
 
 
    // Merge any introduced baffles (from faceZones of faceType 'internal')
    {
        autoPtr<mapPolyMesh> mapPtr = meshRefiner_.mergeZoneBaffles
        (
            true,   // internal zones
            false   // baffle zones
        );
 
        if (mapPtr.valid())
        {
            if (debug & meshRefinement::MESH)
            {
                const_cast<Time&>(mesh.time())++;
                Info<< "Writing baffle-merged mesh to time "
                    << meshRefiner_.timeName() << endl;
                meshRefiner_.write
                (
                    meshRefinement::debugType(debug),
                    meshRefinement::writeType
                    (
                        meshRefinement::writeLevel()
                      | meshRefinement::WRITEMESH
                    ),
                    meshRefiner_.timeName()
                );
            }
        }
    }
 
    // Repatch faces according to nearest. Do not repatch baffle faces.
    {
        labelList duplicateFace(getInternalOrBaffleDuplicateFace());
 
        repatchToSurface(snapParams, adaptPatchIDs, duplicateFace);
    }
 
    if (mergePatchFaces)
    {
        labelList duplicateFace(getInternalOrBaffleDuplicateFace());
 
        // Repatching might have caused faces to be on same patch and hence
        // mergeable so try again to merge coplanar faces. Do not merge baffle
        // faces to ensure they both stay the same.
        label nChanged = meshRefiner_.mergePatchFacesUndo
        (
            featureCos,     // minCos
            featureCos,     // concaveCos
            meshRefiner_.meshedPatches(),
            motionDict,
            duplicateFace   // faces not to merge
        );
 
        nChanged += meshRefiner_.mergeEdgesUndo(featureCos, motionDict);
 
        if (nChanged > 0 && debug & meshRefinement::MESH)
        {
            const_cast<Time&>(mesh.time())++;
            Info<< "Writing patchFace merged mesh to time "
                << meshRefiner_.timeName() << endl;
            meshRefiner_.write
            (
                meshRefinement::debugType(debug),
                meshRefinement::writeType
                (
                    meshRefinement::writeLevel()
                  | meshRefinement::WRITEMESH
                ),
                meshRefiner_.timeName()
            );
        }
    }
 
    if (debug & meshRefinement::MESH)
    {
        const_cast<Time&>(mesh.time())++;
    }
}
 
 
// ************************************************************************* //