From 2c8fa3b35f18be5ff3bb2a491c89229f968cc5d4 Mon Sep 17 00:00:00 2001
From: mattijs <m.janssens@opencfd.co.uk>
Date: Wed, 9 Jul 2008 00:19:33 +0100
Subject: [PATCH] typo

---
 .../autoHexMeshDriverLayers.C                 | 2961 -----------------
 .../autoHexMeshDriverShrink.C                 | 1151 -------
 .../autoHexMeshDriver/autoHexMeshDriverSnap.C | 1245 -------
 .../autoHexMeshDriverTemplates.C              |   75 -
 .../fvMeshDistribute/fvMeshDistribute.H       |    2 +-
 5 files changed, 1 insertion(+), 5433 deletions(-)
 delete mode 100644 src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverLayers.C
 delete mode 100644 src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverShrink.C
 delete mode 100644 src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverSnap.C
 delete mode 100644 src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverTemplates.C

diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverLayers.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverLayers.C
deleted file mode 100644
index 3ebf41390c..0000000000
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverLayers.C
+++ /dev/null
@@ -1,2961 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Description
-    All to do with adding cell layers
-
-\*----------------------------------------------------------------------------*/
-
-#include "autoHexMeshDriver.H"
-#include "fvMesh.H"
-#include "Time.H"
-#include "combineFaces.H"
-#include "removePoints.H"
-#include "pointFields.H"
-#include "motionSmoother.H"
-#include "mathematicalConstants.H"
-#include "pointSet.H"
-#include "faceSet.H"
-#include "cellSet.H"
-#include "polyTopoChange.H"
-#include "mapPolyMesh.H"
-#include "addPatchCellLayer.H"
-#include "mapDistributePolyMesh.H"
-#include "OFstream.H"
-#include "layerParameters.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-Foam::label Foam::autoHexMeshDriver::mergePatchFacesUndo
-(
-    const scalar minCos,
-    const scalar concaveCos,
-    const dictionary& motionDict
-)
-{
-    // Patch face merging engine
-    combineFaces faceCombiner(mesh_, true);
-
-    // Pick up all candidate cells on boundary
-    labelHashSet boundaryCells(mesh_.nFaces()-mesh_.nInternalFaces());
-
-    {
-        labelList patchIDs(meshRefinement::addedPatches(globalToPatch_));
-
-        const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-
-        forAll(patchIDs, i)
-        {
-            label patchI = patchIDs[i];
-
-            const polyPatch& patch = patches[patchI];
-
-            if (!patch.coupled())
-            {
-                forAll(patch, i)
-                {
-                    boundaryCells.insert(mesh_.faceOwner()[patch.start()+i]);
-                }
-            }
-        }
-    }
-
-    // Get all sets of faces that can be merged
-    labelListList allFaceSets
-    (
-        faceCombiner.getMergeSets
-        (
-            minCos,
-            concaveCos,
-            boundaryCells
-        )
-    );
-
-    label nFaceSets = returnReduce(allFaceSets.size(), sumOp<label>());
-
-    Info<< "Merging " << nFaceSets << " sets of faces." << nl << endl;
-
-    if (nFaceSets > 0)
-    {
-        if (debug_)
-        {
-            faceSet allSets(mesh_, "allFaceSets", allFaceSets.size());
-            forAll(allFaceSets, setI)
-            {
-                forAll(allFaceSets[setI], i)
-                {
-                    allSets.insert(allFaceSets[setI][i]);
-                }
-            }
-            Pout<< "Writing all faces to be merged to set "
-                << allSets.objectPath() << endl;
-            allSets.write();
-        }
-
-
-        // Topology changes container
-        polyTopoChange meshMod(mesh_);
-
-        // Merge all faces of a set into the first face of the set.
-        faceCombiner.setRefinement(allFaceSets, meshMod);
-
-        // Experimental: store data for all the points that have been deleted
-        meshRefinerPtr_().storeData
-        (
-            faceCombiner.savedPointLabels(),    // points to store
-            labelList(0),                       // faces to store
-            labelList(0)                        // cells to store
-        );
-
-        // Change the mesh (no inflation)
-        autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
-
-        // Update fields
-        mesh_.updateMesh(map);
-
-        // Move mesh (since morphing does not do this)
-        if (map().hasMotionPoints())
-        {
-            mesh_.movePoints(map().preMotionPoints());
-        }
-
-        faceCombiner.updateMesh(map);
-
-        meshRefinerPtr_().updateMesh(map, labelList(0));
-
-
-
-        for (label iteration = 0; iteration < 100; iteration++)
-        {
-            Info<< nl
-                << "Undo iteration " << iteration << nl
-                << "----------------" << endl;
-
-
-            // Check mesh for errors
-            // ~~~~~~~~~~~~~~~~~~~~~
-
-            faceSet errorFaces
-            (
-                mesh_,
-                "errorFaces",
-                mesh_.nFaces()-mesh_.nInternalFaces()
-            );
-            bool hasErrors = motionSmoother::checkMesh
-            (
-                false,  // report
-                mesh_,
-                motionDict,
-                errorFaces
-            );
-            //if (checkEdgeConnectivity)
-            //{
-            //    Info<< "Checking edge-face connectivity (duplicate faces"
-            //        << " or non-consecutive shared vertices)" << endl;
-            //
-            //    label nOldSize = errorFaces.size();
-            //
-            //    hasErrors =
-            //        mesh_.checkFaceFaces
-            //        (
-            //            false,
-            //            &errorFaces
-            //        )
-            //     || hasErrors;
-            //
-            //    Info<< "Detected additional "
-            //        << returnReduce(errorFaces.size()-nOldSize, sumOp<label>())
-            //        << " faces with illegal face-face connectivity"
-            //        << endl;
-            //}
-
-            if (!hasErrors)
-            {
-                break;
-            }
-
-
-            if (debug_)
-            {
-                Pout<< "Writing all faces in error to faceSet "
-                    << errorFaces.objectPath() << nl << endl;
-                errorFaces.write();
-            }
-
-
-            // Check any master cells for using any of the error faces
-            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-            DynamicList<label> mastersToRestore(allFaceSets.size());
-
-            forAll(allFaceSets, setI)
-            {
-                label masterFaceI = faceCombiner.masterFace()[setI];
-
-                if (masterFaceI != -1)
-                {
-                    label masterCellII = mesh_.faceOwner()[masterFaceI];
-
-                    const cell& cFaces = mesh_.cells()[masterCellII];
-
-                    forAll(cFaces, i)
-                    {
-                        if (errorFaces.found(cFaces[i]))
-                        {
-                            mastersToRestore.append(masterFaceI);
-                            break;
-                        }
-                    }
-                }
-            }
-            mastersToRestore.shrink();
-
-            label nRestore = returnReduce
-            (
-                mastersToRestore.size(),
-                sumOp<label>()
-            );
-
-            Info<< "Masters that need to be restored:"
-                << nRestore << endl;
-
-            if (debug_)
-            {
-                faceSet restoreSet(mesh_, "mastersToRestore", mastersToRestore);
-                Pout<< "Writing all " << mastersToRestore.size()
-                    << " masterfaces to be restored to set "
-                    << restoreSet.objectPath() << endl;
-                restoreSet.write();
-            }
-
-
-            if (nRestore == 0)
-            {
-                break;
-            }
-
-
-            // Undo
-            // ~~~~
-
-            // Topology changes container
-            polyTopoChange meshMod(mesh_);
-
-            // Merge all faces of a set into the first face of the set.
-            // Experimental:mark all points/faces/cells that have been restored.
-            Map<label> restoredPoints(0);
-            Map<label> restoredFaces(0);
-            Map<label> restoredCells(0);
-
-            faceCombiner.setUnrefinement
-            (
-                mastersToRestore,
-                meshMod,
-                restoredPoints,
-                restoredFaces,
-                restoredCells
-            );
-
-            // Change the mesh (no inflation)
-            autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
-
-            // Update fields
-            mesh_.updateMesh(map);
-
-            // Move mesh (since morphing does not do this)
-            if (map().hasMotionPoints())
-            {
-                mesh_.movePoints(map().preMotionPoints());
-            }
-
-            faceCombiner.updateMesh(map);
-
-            // Renumber restore maps
-            inplaceMapKey(map().reversePointMap(), restoredPoints);
-            inplaceMapKey(map().reverseFaceMap(), restoredFaces);
-            inplaceMapKey(map().reverseCellMap(), restoredCells);
-
-            // Experimental:restore all points/face/cells in maps
-            meshRefinerPtr_().updateMesh
-            (
-                map,
-                labelList(0),       // changedFaces
-                restoredPoints,
-                restoredFaces,
-                restoredCells
-            );
-
-            Info<< endl;
-        }
-
-        if (debug_)
-        {
-            Pout<< "Writing merged-faces mesh to time "
-                << mesh_.time().timeName() << nl << endl;
-            mesh_.write();
-        }
-    }
-    else
-    {
-        Info<< "No faces merged ..." << endl;
-    }
-
-    return nFaceSets;
-}
-
-
-// Remove points. pointCanBeDeleted is parallel synchronised.
-Foam::autoPtr<Foam::mapPolyMesh> Foam::autoHexMeshDriver::doRemovePoints
-(
-    removePoints& pointRemover,
-    const boolList& pointCanBeDeleted
-)
-{
-    // Topology changes container
-    polyTopoChange meshMod(mesh_);
-
-    pointRemover.setRefinement(pointCanBeDeleted, meshMod);
-
-    // Change the mesh (no inflation)
-    autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
-
-    // Update fields
-    mesh_.updateMesh(map);
-
-    // Move mesh (since morphing does not do this)
-    if (map().hasMotionPoints())
-    {
-        mesh_.movePoints(map().preMotionPoints());
-    }
-
-    pointRemover.updateMesh(map);
-    meshRefinerPtr_().updateMesh(map, labelList(0));
-
-    return map;
-}
-
-
-// Restore faces (which contain removed points)
-Foam::autoPtr<Foam::mapPolyMesh> Foam::autoHexMeshDriver::doRestorePoints
-(
-    removePoints& pointRemover,
-    const labelList& facesToRestore
-)
-{
-    // Topology changes container
-    polyTopoChange meshMod(mesh_);
-
-    // Determine sets of points and faces to restore
-    labelList localFaces, localPoints;
-    pointRemover.getUnrefimentSet
-    (
-        facesToRestore,
-        localFaces,
-        localPoints
-    );
-
-    // Undo the changes on the faces that are in error.
-    pointRemover.setUnrefinement
-    (
-        localFaces,
-        localPoints,
-        meshMod
-    );
-
-    // Change the mesh (no inflation)
-    autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
-
-    // Update fields
-    mesh_.updateMesh(map);
-
-    // Move mesh (since morphing does not do this)
-    if (map().hasMotionPoints())
-    {
-        mesh_.movePoints(map().preMotionPoints());
-    }
-
-    pointRemover.updateMesh(map);
-    meshRefinerPtr_().updateMesh(map, labelList(0));
-
-    return map;
-}
-
-
-// Collect all faces that are both in candidateFaces and in the set.
-// If coupled face also collects the coupled face.
-Foam::labelList Foam::autoHexMeshDriver::collectFaces
-(
-    const labelList& candidateFaces,
-    const labelHashSet& set
-) const
-{
-    // Has face been selected?
-    boolList selected(mesh_.nFaces(), false);
-
-    forAll(candidateFaces, i)
-    {
-        label faceI = candidateFaces[i];
-
-        if (set.found(faceI))
-        {
-            selected[faceI] = true;
-        }
-    }
-    syncTools::syncFaceList
-    (
-        mesh_,
-        selected,
-        orEqOp<bool>(),     // combine operator
-        false               // separation
-    );
-
-    labelList selectedFaces(findIndices(selected, true));
-
-    return selectedFaces;
-}
-
-
-// Pick up faces of cells of faces in set.
-Foam::labelList Foam::autoHexMeshDriver::growFaceCellFace
-(
-    const labelHashSet& set
-) const
-{
-    boolList selected(mesh_.nFaces(), false);
-
-    forAllConstIter(faceSet, set, iter)
-    {
-        label faceI = iter.key();
-
-        label own = mesh_.faceOwner()[faceI];
-
-        const cell& ownFaces = mesh_.cells()[own];
-        forAll(ownFaces, ownFaceI)
-        {
-            selected[ownFaces[ownFaceI]] = true;
-        }
-
-        if (mesh_.isInternalFace(faceI))
-        {
-            label nbr = mesh_.faceNeighbour()[faceI];
-
-            const cell& nbrFaces = mesh_.cells()[nbr];
-            forAll(nbrFaces, nbrFaceI)
-            {
-                selected[nbrFaces[nbrFaceI]] = true;
-            }
-        }
-    }
-    syncTools::syncFaceList
-    (
-        mesh_,
-        selected,
-        orEqOp<bool>(),     // combine operator
-        false               // separation
-    );
-    return findIndices(selected, true);
-}
-
-
-// Remove points not used by any face or points used by only two faces where
-// the edges are in line
-Foam::label Foam::autoHexMeshDriver::mergeEdgesUndo
-(
-    const scalar minCos,
-    const dictionary& motionDict
-)
-{
-    Info<< nl
-        << "Merging all points on surface that" << nl
-        << "- are used by only two boundary faces and" << nl
-        << "- make an angle with a cosine of more than " << minCos
-        << "." << nl << endl;
-
-    // Point removal analysis engine with undo
-    removePoints pointRemover(mesh_, true);
-
-    // Count usage of points
-    boolList pointCanBeDeleted;
-    label nRemove = pointRemover.countPointUsage(minCos, pointCanBeDeleted);
-
-    if (nRemove > 0)
-    {
-        Info<< "Removing " << nRemove
-            << " straight edge points ..." << nl << endl;
-
-        // Remove points
-        // ~~~~~~~~~~~~~
-
-        doRemovePoints(pointRemover, pointCanBeDeleted);
-
-
-        for (label iteration = 0; iteration < 100; iteration++)
-        {
-            Info<< nl
-                << "Undo iteration " << iteration << nl
-                << "----------------" << endl;
-
-
-            // Check mesh for errors
-            // ~~~~~~~~~~~~~~~~~~~~~
-
-            faceSet errorFaces
-            (
-                mesh_,
-                "errorFaces",
-                mesh_.nFaces()-mesh_.nInternalFaces()
-            );
-            bool hasErrors = motionSmoother::checkMesh
-            (
-                false,  // report
-                mesh_,
-                motionDict,
-                errorFaces
-            );
-            //if (checkEdgeConnectivity)
-            //{
-            //    Info<< "Checking edge-face connectivity (duplicate faces"
-            //        << " or non-consecutive shared vertices)" << endl;
-            //
-            //    label nOldSize = errorFaces.size();
-            //
-            //    hasErrors =
-            //        mesh_.checkFaceFaces
-            //        (
-            //            false,
-            //            &errorFaces
-            //        )
-            //     || hasErrors;
-            //
-            //    Info<< "Detected additional "
-            //        << returnReduce(errorFaces.size()-nOldSize, sumOp<label>())
-            //        << " faces with illegal face-face connectivity"
-            //        << endl;
-            //}
-
-            if (!hasErrors)
-            {
-                break;
-            }
-
-            if (debug_)
-            {
-                Pout<< "**Writing all faces in error to faceSet "
-                    << errorFaces.objectPath() << nl << endl;
-                errorFaces.write();
-            }
-
-            labelList masterErrorFaces
-            (
-                collectFaces
-                (
-                    pointRemover.savedFaceLabels(),
-                    errorFaces
-                )
-            );
-
-            label n = returnReduce(masterErrorFaces.size(), sumOp<label>());
-
-            Info<< "Detected " << n
-                << " error faces on boundaries that have been merged."
-                << " These will be restored to their original faces." << nl
-                << endl;
-
-            if (n == 0)
-            {
-                if (hasErrors)
-                {
-                    Info<< "Detected "
-                        << returnReduce(errorFaces.size(), sumOp<label>())
-                        << " error faces in mesh."
-                        << " Restoring neighbours of faces in error." << nl
-                        << endl;
-
-                    labelList expandedErrorFaces
-                    (
-                        growFaceCellFace
-                        (
-                            errorFaces
-                        )
-                    );
-
-                    doRestorePoints(pointRemover, expandedErrorFaces);
-                }
-
-                break;
-            }
-
-            doRestorePoints(pointRemover, masterErrorFaces);
-        }
-
-        if (debug_)
-        {
-            Pout<< "Writing merged-edges mesh to time "
-                << mesh_.time().timeName() << nl << endl;
-            mesh_.write();
-        }
-    }
-    else
-    {
-        Info<< "No straight edges simplified and no points removed ..." << endl;
-    }
-
-    return nRemove;
-}
-
-
-// For debugging: Dump displacement to .obj files
-void Foam::autoHexMeshDriver::dumpDisplacement
-(
-    const fileName& prefix,
-    const indirectPrimitivePatch& pp,
-    const vectorField& patchDisp,
-    const List<extrudeMode>& extrudeStatus
-)
-{
-    OFstream dispStr(prefix + "_disp.obj");
-    Info<< "Writing all displacements to " << dispStr.name() << nl << endl;
-
-    label vertI = 0;
-
-    forAll(patchDisp, patchPointI)
-    {
-        const point& pt = pp.localPoints()[patchPointI];
-
-        meshTools::writeOBJ(dispStr, pt); vertI++;
-        meshTools::writeOBJ(dispStr, pt + patchDisp[patchPointI]); vertI++;
-
-        dispStr << "l " << vertI-1 << ' ' << vertI << nl;
-    }
-
-
-    OFstream illStr(prefix + "_illegal.obj");
-    Info<< "Writing invalid displacements to " << illStr.name() << nl << endl;
-
-    vertI = 0;
-
-    forAll(patchDisp, patchPointI)
-    {
-        if (extrudeStatus[patchPointI] != EXTRUDE)
-        {
-            const point& pt = pp.localPoints()[patchPointI];
-
-            meshTools::writeOBJ(illStr, pt); vertI++;
-            meshTools::writeOBJ(illStr, pt + patchDisp[patchPointI]); vertI++;
-
-            illStr << "l " << vertI-1 << ' ' << vertI << nl;
-        }
-    }
-}
-
-
-// Check that primitivePatch is not multiply connected. Collect non-manifold
-// points in pointSet.
-void Foam::autoHexMeshDriver::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::autoHexMeshDriver::checkMeshManifold() const
-{
-    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.write();
-    }
-    Info<< endl;
-}
-
-
-
-// Unset extrusion on point. Returns true if anything unset.
-bool Foam::autoHexMeshDriver::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::autoHexMeshDriver::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::autoHexMeshDriver::handleNonManifolds
-(
-    const indirectPrimitivePatch& pp,
-    const labelList& meshEdges,
-    pointField& patchDisp,
-    labelList& patchNLayers,
-    List<extrudeMode>& extrudeStatus
-) const
-{
-    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);
-
-    label nNonManif = returnReduce(nonManifoldPoints.size(), sumOp<label>());
-
-    Info<< "Outside of local patch is multiply connected across edges or"
-        << " points at " << nNonManif << " points." << endl;
-
-
-    const labelList& meshPoints = pp.meshPoints();
-
-
-    // 2. Parallel check
-    // For now disable extrusion at any shared edges.
-    const labelHashSet sharedEdgeSet(mesh_.globalData().sharedEdgeLabels());
-
-    forAll(pp.edges(), edgeI)
-    {
-        if (sharedEdgeSet.found(meshEdges[edgeI]))
-        {
-            const edge& e = mesh_.edges()[meshEdges[edgeI]];
-
-            Pout<< "Disabling extrusion at edge "
-                << mesh_.points()[e[0]] << mesh_.points()[e[1]]
-                << " since it is non-manifold across coupled patches."
-                << endl;
-
-            nonManifoldPoints.insert(e[0]);
-            nonManifoldPoints.insert(e[1]);
-        }
-    }
-
-    // 3b. extrusion can produce multiple faces between 2 cells
-    // across processor boundary
-    // This occurs when a coupled face shares more than 1 edge with a
-    // non-coupled boundary face.
-    // This is now correctly handled by addPatchCellLayer in that it
-    // extrudes a single face from the stringed up edges.
-
-
-    nNonManif = returnReduce(nonManifoldPoints.size(), sumOp<label>());
-
-    if (nNonManif > 0)
-    {
-        Info<< "Outside of patches is multiply connected across edges or"
-            << " points at " << nNonManif << " points." << nl
-            << "Writing " << nNonManif
-            << " points where this happens to pointSet "
-            << nonManifoldPoints.name() << nl
-            << "and setting layer thickness to zero on these points."
-            << endl;
-
-        nonManifoldPoints.write();
-
-        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::autoHexMeshDriver::handleFeatureAngle
-(
-    const indirectPrimitivePatch& pp,
-    const labelList& meshEdges,
-    const scalar minCos,
-    pointField& patchDisp,
-    labelList& patchNLayers,
-    List<extrudeMode>& extrudeStatus
-) const
-{
-    Info<< nl << "Handling feature edges ..." << endl;
-
-    if (minCos < 1-SMALL)
-    {
-        // Normal component of normals of connected faces.
-        vectorField edgeNormal(mesh_.nEdges(), wallPoint::greatPoint);
-
-        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(),
-            wallPoint::greatPoint,  // null value
-            false                   // no separation
-        );
-
-        label vertI = 0;
-        autoPtr<OFstream> str;
-        if (debug_)
-        {
-            str.reset(new OFstream(mesh_.time().path()/"featureEdges.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 != wallPoint::greatPoint)
-            {
-                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())
-                    {
-                        meshTools::writeOBJ(str(), pp.localPoints()[e[0]]);
-                        vertI++;
-                        meshTools::writeOBJ(str(), pp.localPoints()[e[1]]);
-                        vertI++;
-                        str()<< "l " << vertI-1 << ' ' << vertI << nl;
-                    }
-                }
-            }
-        }
-
-        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::autoHexMeshDriver::handleWarpedFaces
-(
-    const indirectPrimitivePatch& pp,
-    const scalar faceRatio,
-    const scalar edge0Len,
-    const labelList& cellLevel,
-    pointField& patchDisp,
-    labelList& patchNLayers,
-    List<extrudeMode>& extrudeStatus
-) const
-{
-    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::autoHexMeshDriver::handleMultiplePatchFaces
-//(
-//    const indirectPrimitivePatch& pp,
-//    pointField& patchDisp,
-//    labelList& patchNLayers,
-//    List<extrudeMode>& extrudeStatus
-//) const
-//{
-//    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)
-//    {
-//        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;
-//}
-
-
-// No extrusion on faces with differing number of layers for points
-void Foam::autoHexMeshDriver::setNumLayers
-(
-    const labelList& patchToNLayers,
-    const labelList& patchIDs,
-    const indirectPrimitivePatch& pp,
-    pointField& patchDisp,
-    labelList& patchNLayers,
-    List<extrudeMode>& extrudeStatus
-) const
-{
-    Info<< nl << "Handling points with inconsistent layer specification ..."
-        << endl;
-
-    // Get for every point (really only nessecary 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
-        false               // no separation
-    );
-    syncTools::syncPointList
-    (
-        mesh_,
-        pp.meshPoints(),
-        minLayers,
-        minEqOp<label>(),
-        labelMax,           // null value
-        false               // no separation
-    );
-
-    // 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];
-        }
-    }
-
-    //reduce(nConflicts, sumOp<label>());
-    //
-    //Info<< "Set displacement to zero for " << nConflicts
-    //    << " points due to points being on multiple regions"
-    //    << " with inconsistent nLayers specification." << endl;
-}
-
-
-// Grow no-extrusion layer
-void Foam::autoHexMeshDriver::growNoExtrusion
-(
-    const indirectPrimitivePatch& pp,
-    pointField& patchDisp,
-    labelList& patchNLayers,
-    List<extrudeMode>& extrudeStatus
-)
-{
-    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]] == NOEXTRUDE
-                 && grownExtrudeStatus[f[fp]] != NOEXTRUDE
-                )
-                {
-                    grownExtrudeStatus[f[fp]] = NOEXTRUDE;
-                    nGrown++;
-                }
-            }
-        }
-    }
-
-    extrudeStatus.transfer(grownExtrudeStatus);
-
-    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::autoHexMeshDriver::calculateLayerThickness
-(
-    const indirectPrimitivePatch& pp,
-    const labelList& patchIDs,
-    const scalarField& patchExpansionRatio,
-    const scalarField& patchFinalLayerRatio,
-    const scalarField& patchRelMinThickness,
-    const labelList& cellLevel,
-    const labelList& patchNLayers,
-    const scalar edge0Len,
-
-    scalarField& thickness,
-    scalarField& minThickness,
-    scalarField& expansionRatio
-) const
-{
-    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-
-    if (min(patchRelMinThickness) < 0 || max(patchRelMinThickness) > 2)
-    {
-        FatalErrorIn("calculateLayerThickness(..)")
-            << "Thickness should be factor of local undistorted cell size."
-            << " Valid values are [0..2]." << nl
-            << " minThickness:" << patchRelMinThickness
-            << exit(FatalError);
-    }
-
-
-    // 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
-            false               // no separation
-        );
-    }
-
-
-    // Rework patch-wise layer parameters into minimum per point
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    expansionRatio.setSize(pp.nPoints());
-    expansionRatio = GREAT;
-    scalarField finalLayerRatio(pp.nPoints(), GREAT);
-    scalarField relMinThickness(pp.nPoints(), GREAT);
-
-    {
-        forAll(patchIDs, i)
-        {
-            label patchI = patchIDs[i];
-
-            const labelList& meshPoints = patches[patchI].meshPoints();
-
-            forAll(meshPoints, patchPointI)
-            {
-                label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];
-
-                expansionRatio[ppPointI] = min
-                (
-                    expansionRatio[ppPointI],
-                    patchExpansionRatio[patchI]
-                );
-                finalLayerRatio[ppPointI] = min
-                (
-                    finalLayerRatio[ppPointI],
-                    patchFinalLayerRatio[patchI]
-                );
-                relMinThickness[ppPointI] = min
-                (
-                    relMinThickness[ppPointI],
-                    patchRelMinThickness[patchI]
-                );
-            }
-        }
-
-        syncTools::syncPointList
-        (
-            mesh_,
-            pp.meshPoints(),
-            expansionRatio,
-            minEqOp<scalar>(),
-            GREAT,              // null value
-            false               // no separation
-        );
-        syncTools::syncPointList
-        (
-            mesh_,
-            pp.meshPoints(),
-            finalLayerRatio,
-            minEqOp<scalar>(),
-            GREAT,              // null value
-            false               // no separation
-        );
-        syncTools::syncPointList
-        (
-            mesh_,
-            pp.meshPoints(),
-            relMinThickness,
-            minEqOp<scalar>(),
-            GREAT,              // null value
-            false               // no separation
-        );
-    }
-
-
-
-    // Per mesh point the expansion parameters
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    thickness.setSize(pp.nPoints());
-    minThickness.setSize(pp.nPoints());
-
-    forAll(maxPointLevel, pointI)
-    {
-        // Find undistorted edge size for this level.
-        scalar edgeLen = edge0Len/(1<<maxPointLevel[pointI]);
-
-        // Calculate layer thickness based on expansion ratio
-        // and final layer height
-        if (expansionRatio[pointI] == 1)
-        {
-            thickness[pointI] =
-                finalLayerRatio[pointI]
-              * patchNLayers[pointI]
-              * edgeLen;
-            minThickness[pointI] = relMinThickness[pointI]*edgeLen;
-        }
-        else
-        {
-            scalar invExpansion = 1.0 / expansionRatio[pointI];
-            label nLay = patchNLayers[pointI];
-            thickness[pointI] =
-                finalLayerRatio[pointI]
-              * edgeLen
-              * (1.0 - pow(invExpansion, nLay))
-              / (1.0 - invExpansion);
-            minThickness[pointI] = relMinThickness[pointI]*edgeLen;
-        }
-    }
-
-    Info<< "calculateLayerThickness : min:" << gMin(thickness)
-        << " max:" << gMax(thickness) << endl;
-}
-
-
-// Synchronize displacement among coupled patches.
-void Foam::autoHexMeshDriver::syncPatchDisplacement
-(
-    const motionSmoother& meshMover,
-    const scalarField& minThickness,
-    pointField& patchDisp,
-    labelList& patchNLayers,
-    List<extrudeMode>& extrudeStatus
-) const
-{
-    const labelList& meshPoints = meshMover.patch().meshPoints();
-
-    label nChangedTotal = 0;
-
-    while (true)
-    {
-        label nChanged = 0;
-
-        // Sync displacement (by taking min)
-        syncTools::syncPointList
-        (
-            mesh_,
-            meshPoints,
-            patchDisp,
-            minEqOp<vector>(),
-            wallPoint::greatPoint,          // null value
-            false                           // no separation
-        );
-
-        // 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
-            false               // no separation
-        );
-
-        // Reset if differs
-        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
-            false               // no separation
-        );
-
-        // Reset if differs
-        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::autoHexMeshDriver::getPatchDisplacement
-(
-    const motionSmoother& meshMover,
-    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 indirectPrimitivePatch& pp = meshMover.patch();
-    const vectorField& faceNormals = pp.faceNormals();
-    const labelListList& pointFaces = pp.pointFaces();
-    const pointField& localPoints = pp.localPoints();
-    const labelList& meshPoints = pp.meshPoints();
-
-    // Determine pointNormal
-    // ~~~~~~~~~~~~~~~~~~~~~
-
-    pointField pointNormals(pp.nPoints(), vector::zero);
-    {
-        labelList nPointFaces(pp.nPoints(), 0);
-
-        forAll(faceNormals, faceI)
-        {
-            const face& f = pp.localFaces()[faceI];
-
-            forAll(f, fp)
-            {
-                pointNormals[f[fp]] += faceNormals[faceI];
-                nPointFaces[f[fp]] ++;
-            }
-        }
-
-        syncTools::syncPointList
-        (
-            mesh_,
-            meshPoints,
-            pointNormals,
-            plusEqOp<vector>(),
-            vector::zero,       // null value
-            false               // no separation
-        );
-
-        syncTools::syncPointList
-        (
-            mesh_,
-            meshPoints,
-            nPointFaces,
-            plusEqOp<label>(),
-            0,                  // null value
-            false               // no separation
-        );
-
-        forAll(pointNormals, i)
-        {
-            pointNormals[i] /= nPointFaces[i];
-        }
-    }
-
-
-    // Determine local length scale on patch
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Start off from same thickness everywhere (except where no extrusion)
-    patchDisp = thickness*pointNormals;
-
-    // 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]))
-            {
-                Pout<< "No valid normal for point " << meshPointI
-                    << ' ' << pp.points()[meshPointI]
-                    << "; setting displacement to " << patchDisp[patchPointI]
-                    << endl;
-
-                extrudeStatus[patchPointI] = EXTRUDEREMOVE;
-            }
-        }
-    }
-
-    // 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)
-            {
-                Pout<< "Displacement at illegal point "
-                    << localPoints[patchPointI]
-                    << " set to " << (avg / nPoints - localPoints[patchPointI])
-                    << endl;
-
-                patchDisp[patchPointI] =
-                    avg / nPoints
-                  - localPoints[patchPointI];
-            }
-        }
-    }
-
-    // Make sure displacement is equal on both sides of coupled patches.
-    syncPatchDisplacement
-    (
-        meshMover,
-        minThickness,
-        patchDisp,
-        patchNLayers,
-        extrudeStatus
-    );
-
-    Info<< endl;
-}
-
-
-// Truncates displacement
-// - for all patchFaces in the faceset displacement gets set to zero
-// - all displacement < minThickness gets set to zero
-Foam::label Foam::autoHexMeshDriver::truncateDisplacement
-(
-    const motionSmoother& meshMover,
-    const scalarField& minThickness,
-    const faceSet& illegalPatchFaces,
-    pointField& patchDisp,
-    labelList& patchNLayers,
-    List<extrudeMode>& extrudeStatus
-) const
-{
-    const polyMesh& mesh = meshMover.mesh();
-    const indirectPrimitivePatch& pp = meshMover.patch();
-
-    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
-        (
-            meshMover,
-            minThickness,
-            patchDisp,
-            patchNLayers,
-            extrudeStatus
-        );
-
-        // 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;
-
-
-        // Make sure that a face has consistent number of layers for all
-        // its vertices.
-
-        label nDiffering = 0;
-
-        //forAll(localFaces, i)
-        //{
-        //    const face& localF = localFaces[i];
-        //
-        //    label numLayers = -1;
-        //
-        //    forAll(localF, fp)
-        //    {
-        //        if (patchNLayers[localF[fp]] > 0)
-        //        {
-        //            if (numLayers == -1)
-        //            {
-        //                numLayers = patchNLayers[localF[fp]];
-        //            }
-        //            else if (numLayers != patchNLayers[localF[fp]])
-        //            {
-        //                // Differing number of layers
-        //                if
-        //                (
-        //                    unmarkExtrusion
-        //                    (
-        //                        localF,
-        //                        patchDisp,
-        //                        patchNLayers,
-        //                        extrudeStatus
-        //                    )
-        //                )
-        //                {
-        //                    nDiffering++;
-        //                    nChanged++;
-        //                }
-        //                break;
-        //            }
-        //        }
-        //    }
-        //}
-        //
-        //reduce(nDiffering, sumOp<label>());
-        //
-        //Info<< "truncateDisplacement : Unextruded " << nDiffering
-        //    << " faces due to having differing number of layers." << endl;
-
-        if (nPinched+nDiffering == 0)
-        {
-            break;
-        }
-    }
-
-    return nChanged;
-}
-
-
-// Setup layer information (at points and faces) to modify mesh topology in
-// regions where layer mesh terminates.
-void Foam::autoHexMeshDriver::setupLayerInfoTruncation
-(
-    const motionSmoother& meshMover,
-    const labelList& patchNLayers,
-    const List<extrudeMode>& extrudeStatus,
-    const label nBufferCellsNoExtrude,
-    labelList& nPatchPointLayers,
-    labelList& nPatchFaceLayers
-) const
-{
-    Info<< nl << "Setting up information for layer truncation ..." << endl;
-
-    const indirectPrimitivePatch& pp = meshMover.patch();
-    const polyMesh& mesh = meshMover.mesh();
-
-    if (nBufferCellsNoExtrude < 0)
-    {
-        Info<< nl << "Performing no layer truncation."
-            << " nBufferCellsNoExtrude set to less than 0  ..." << endl;
-
-        // Face layers if any point get extruded
-        forAll(pp.localFaces(), patchFaceI)
-        {
-            const face& f = pp.localFaces()[patchFaceI];
-
-            forAll(f, fp)
-            {
-                if (patchNLayers[f[fp]] > 0)
-                {
-                    nPatchFaceLayers[patchFaceI] = patchNLayers[f[fp]];
-                    break;
-                }
-            }
-        }
-        nPatchPointLayers = patchNLayers;
-    }
-    else
-    {
-        // Determine max point layers per face.
-        labelList maxLevel(pp.size(), 0);
-
-        forAll(pp.localFaces(), patchFaceI)
-        {
-            const face& f = pp.localFaces()[patchFaceI];
-
-            // find patch faces where layer terminates (i.e contains extrude
-            // and noextrude points).
-
-            bool noExtrude = false;
-            label mLevel = 0;
-
-            forAll(f, fp)
-            {
-                if (extrudeStatus[f[fp]] == NOEXTRUDE)
-                {
-                    noExtrude = true;
-                }
-                mLevel = max(mLevel, patchNLayers[f[fp]]);
-            }
-
-            if (mLevel > 0)
-            {
-                // So one of the points is extruded. Check if all are extruded
-                // or is a mix.
-
-                if (noExtrude)
-                {
-                    nPatchFaceLayers[patchFaceI] = 1;
-                    maxLevel[patchFaceI] = mLevel;
-                }
-                else
-                {
-                    maxLevel[patchFaceI] = mLevel;
-                }
-            }
-        }
-
-        // We have the seed faces (faces with nPatchFaceLayers != maxLevel)
-        // Now do a meshwave across the patch where we pick up neighbours
-        // of seed faces.
-        // Note: quite inefficient. Could probably be coded better.
-
-        const labelListList& pointFaces = pp.pointFaces();
-
-        label nLevels = gMax(patchNLayers);
-
-        // flag neighbouring patch faces with number of layers to grow
-        for (label ilevel = 1; ilevel < nLevels; ilevel++)
-        {
-            label nBuffer;
-
-            if (ilevel == 1)
-            {
-                nBuffer = nBufferCellsNoExtrude - 1;
-            }
-            else
-            {
-                nBuffer = nBufferCellsNoExtrude;
-            }
-
-            for (label ibuffer = 0; ibuffer < nBuffer + 1; ibuffer++)
-            {
-                labelList tempCounter(nPatchFaceLayers);
-
-                boolList foundNeighbour(pp.nPoints(), false);
-
-                forAll(pp.meshPoints(), patchPointI)
-                {
-                    forAll(pointFaces[patchPointI], pointFaceI)
-                    {
-                        label faceI = pointFaces[patchPointI][pointFaceI];
-
-                        if
-                        (
-                            nPatchFaceLayers[faceI] != -1
-                         && maxLevel[faceI] > 0
-                        )
-                        {
-                            foundNeighbour[patchPointI] = true;
-                            break;
-                        }
-                    }
-                }
-
-                syncTools::syncPointList
-                (
-                    mesh,
-                    pp.meshPoints(),
-                    foundNeighbour,
-                    orEqOp<bool>(),
-                    false,              // null value
-                    false               // no separation
-                );
-
-                forAll(pp.meshPoints(), patchPointI)
-                {
-                    if (foundNeighbour[patchPointI])
-                    {
-                        forAll(pointFaces[patchPointI], pointFaceI)
-                        {
-                            label faceI = pointFaces[patchPointI][pointFaceI];
-                            if
-                            (
-                                nPatchFaceLayers[faceI] == -1
-                             && maxLevel[faceI] > 0
-                             && ilevel < maxLevel[faceI]
-                            )
-                            {
-                                tempCounter[faceI] = ilevel;
-                            }
-                        }
-                    }
-                }
-                nPatchFaceLayers = tempCounter;
-            }
-        }
-
-        forAll(pp.localFaces(), patchFaceI)
-        {
-            if (nPatchFaceLayers[patchFaceI] == -1)
-            {
-                nPatchFaceLayers[patchFaceI] = maxLevel[patchFaceI];
-            }
-        }
-
-        forAll(pp.meshPoints(), patchPointI)
-        {
-            if (extrudeStatus[patchPointI] != NOEXTRUDE)
-            {
-                forAll(pointFaces[patchPointI], pointFaceI)
-                {
-                    label face = pointFaces[patchPointI][pointFaceI];
-                    nPatchPointLayers[patchPointI] = max
-                    (
-                        nPatchPointLayers[patchPointI],
-                        nPatchFaceLayers[face]
-                    );
-                }
-            }
-            else
-            {
-                nPatchPointLayers[patchPointI] = 0;
-            }
-        }
-        syncTools::syncPointList
-        (
-            mesh,
-            pp.meshPoints(),
-            nPatchPointLayers,
-            maxEqOp<label>(),
-            0,                  // null value
-            false               // no separation
-        );
-    }
-}
-
-
-// Does any of the cells use a face from faces?
-bool Foam::autoHexMeshDriver::cellsUseFace
-(
-    const polyMesh& mesh,
-    const labelList& cellLabels,
-    const labelHashSet& faces
-)
-{
-    forAll(cellLabels, i)
-    {
-        const cell& cFaces = mesh.cells()[cellLabels[i]];
-
-        forAll(cFaces, cFaceI)
-        {
-            if (faces.found(cFaces[cFaceI]))
-            {
-                return true;
-            }
-        }
-    }
-    return false;
-}
-
-
-// Checks the newly added cells and locally unmarks points so they
-// will not get extruded next time round. Returns global number of unmarked
-// points (0 if all was fine)
-Foam::label Foam::autoHexMeshDriver::checkAndUnmark
-(
-    const addPatchCellLayer& addLayer,
-    const dictionary& motionDict,
-    const indirectPrimitivePatch& pp,
-    const fvMesh& newMesh,
-
-    pointField& patchDisp,
-    labelList& patchNLayers,
-    List<extrudeMode>& extrudeStatus
-)
-{
-    // Check the resulting mesh for errors
-    Info<< nl << "Checking mesh with layer ..." << endl;
-    faceSet wrongFaces(newMesh, "wrongFaces", newMesh.nFaces()/1000);
-    motionSmoother::checkMesh(false, newMesh, motionDict, wrongFaces);
-    Info<< "Detected " << returnReduce(wrongFaces.size(), sumOp<label>())
-        << " illegal faces"
-        << " (concave, zero area or negative cell pyramid volume)"
-        << endl;
-
-    // Undo local extrusion if
-    // - any of the added cells in error
-
-    label nChanged = 0;
-
-    // Get all cells in the layer.
-    labelListList addedCells
-    (
-        addPatchCellLayer::addedCells
-        (
-            newMesh,
-            addLayer.layerFaces()
-        )
-    );
-
-    // Check if any of the faces in error uses any face of an added cell
-    forAll(addedCells, oldPatchFaceI)
-    {
-        // Get the cells (in newMesh labels) per old patch face (in mesh
-        // labels)
-        const labelList& fCells = addedCells[oldPatchFaceI];
-
-        if (cellsUseFace(newMesh, fCells, wrongFaces))
-        {
-            // Unmark points on old mesh
-            if
-            (
-                unmarkExtrusion
-                (
-                    pp.localFaces()[oldPatchFaceI],
-                    patchDisp,
-                    patchNLayers,
-                    extrudeStatus
-                )
-            )
-            {
-                nChanged++;
-            }
-        }
-    }
-
-    return returnReduce(nChanged, sumOp<label>());
-}
-
-
-//- Count global number of extruded faces
-Foam::label Foam::autoHexMeshDriver::countExtrusion
-(
-    const indirectPrimitivePatch& pp,
-    const List<extrudeMode>& extrudeStatus
-)
-{
-    // Count number of extruded patch faces
-    label nExtruded = 0;
-    {
-        const faceList& localFaces = pp.localFaces();
-
-        forAll(localFaces, i)
-        {
-            const face& localFace = localFaces[i];
-
-            forAll(localFace, fp)
-            {
-                if (extrudeStatus[localFace[fp]] != NOEXTRUDE)
-                {
-                    nExtruded++;
-                    break;
-                }
-            }
-        }
-    }
-
-    return returnReduce(nExtruded, sumOp<label>());
-}
-
-
-// Collect layer faces and layer cells into bools for ease of handling
-void Foam::autoHexMeshDriver::getLayerCellsFaces
-(
-    const polyMesh& mesh,
-    const addPatchCellLayer& addLayer,
-    boolList& flaggedCells,
-    boolList& flaggedFaces
-)
-{
-    flaggedCells.setSize(mesh.nCells());
-    flaggedCells = false;
-    flaggedFaces.setSize(mesh.nFaces());
-    flaggedFaces = false;
-
-    // Mark all faces in the layer
-    const labelListList& layerFaces = addLayer.layerFaces();
-
-    // Mark all cells in the layer.
-    labelListList addedCells(addPatchCellLayer::addedCells(mesh, layerFaces));
-
-    forAll(addedCells, oldPatchFaceI)
-    {
-        const labelList& added = addedCells[oldPatchFaceI];
-
-        forAll(added, i)
-        {
-            flaggedCells[added[i]] = true;
-        }
-    }
-
-    forAll(layerFaces, oldPatchFaceI)
-    {
-        const labelList& layer = layerFaces[oldPatchFaceI];
-
-        if (layer.size() > 0)
-        {
-            for (label i = 1; i < layer.size()-1; i++)
-            {
-                flaggedFaces[layer[i]] = true;
-            }
-        }
-    }
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::autoHexMeshDriver::mergePatchFacesUndo
-(
-    const layerParameters& layerParams,
-    const dictionary& motionDict
-)
-{
-    scalar minCos = Foam::cos
-    (
-        layerParams.featureAngle()
-      * mathematicalConstant::pi/180.0
-    );
-
-    scalar concaveCos = Foam::cos
-    (
-        layerParams.concaveAngle()
-      * mathematicalConstant::pi/180.0
-    );
-
-    Info<< nl
-        << "Merging all faces of a cell" << nl
-        << "---------------------------" << nl
-        << "    - which are on the same patch" << nl
-        << "    - which make an angle < " << layerParams.featureAngle()
-        << " degrees"
-        << nl
-        << "      (cos:" << minCos << ')' << nl
-        << "    - as long as the resulting face doesn't become concave"
-        << " by more than "
-        << layerParams.concaveAngle()
-        << " degrees (0=straight, 180=fully concave)" << nl
-        << endl;
-
-    label nChanged = mergePatchFacesUndo(minCos, concaveCos, motionDict);
-
-    nChanged += mergeEdgesUndo(minCos, motionDict);
-}
-
-
-void Foam::autoHexMeshDriver::addLayers
-(
-    const layerParameters& layerParams,
-    const dictionary& motionDict,
-    const label nAllowableErrors,
-    motionSmoother& meshMover
-)
-{
-    const indirectPrimitivePatch& pp = meshMover.patch();
-    const labelList& meshPoints = pp.meshPoints();
-
-    // Precalculate mesh edge labels for patch edges
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    labelList meshEdges(pp.nEdges());
-    forAll(pp.edges(), edgeI)
-    {
-        const edge& ppEdge = pp.edges()[edgeI];
-        label v0 = meshPoints[ppEdge[0]];
-        label v1 = meshPoints[ppEdge[1]];
-        meshEdges[edgeI] = meshTools::findEdge
-        (
-            mesh_.edges(),
-            mesh_.pointEdges()[v0],
-            v0,
-            v1
-        );
-    }
-
-    // Displacement for all pp.localPoints.
-    vectorField patchDisp(pp.nPoints(), vector::one);
-
-    // Number of layers for all pp.localPoints. Note: only valid if
-    // extrudeStatus = EXTRUDE.
-    labelList patchNLayers(pp.nPoints(), 0);
-
-    // Whether to add edge for all pp.localPoints.
-    List<extrudeMode> extrudeStatus(pp.nPoints(), EXTRUDE);
-
-
-    // Get number of layer per point from number of layers per patch
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    setNumLayers
-    (
-        layerParams.numLayers(),    // per patch the num layers
-        meshMover.adaptPatchIDs(),  // patches that are being moved
-        pp,                         // indirectpatch for all faces moving
-
-        patchDisp,
-        patchNLayers,
-        extrudeStatus
-    );
-
-    // Disable extrusion on non-manifold points
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    handleNonManifolds
-    (
-        pp,
-        meshEdges,
-
-        patchDisp,
-        patchNLayers,
-        extrudeStatus
-    );
-
-    // Disable extrusion on feature angles
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    handleFeatureAngle
-    (
-        pp,
-        meshEdges,
-        layerParams.featureAngle()*mathematicalConstant::pi/180.0,
-
-        patchDisp,
-        patchNLayers,
-        extrudeStatus
-    );
-
-    // Disable extrusion on warped faces
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Undistorted edge length
-    const scalar edge0Len = meshRefinerPtr_().meshCutter().level0EdgeLength();
-    const labelList& cellLevel = meshRefinerPtr_().meshCutter().cellLevel();
-
-    handleWarpedFaces
-    (
-        pp,
-        layerParams.maxFaceThicknessRatio(),
-        edge0Len,
-        cellLevel,
-
-        patchDisp,
-        patchNLayers,
-        extrudeStatus
-    );
-
-    //// Disable extrusion on cells with multiple patch faces
-    //// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    //
-    //handleMultiplePatchFaces
-    //(
-    //    pp,
-    //
-    //    patchDisp,
-    //    patchNLayers,
-    //    extrudeStatus
-    //);
-
-
-    // Grow out region of non-extrusion
-    for (label i = 0; i < layerParams.nGrow(); i++)
-    {
-        growNoExtrusion
-        (
-            pp,
-            patchDisp,
-            patchNLayers,
-            extrudeStatus
-        );
-    }
-
-    // Determine (wanted) point-wise layer thickness and expansion ratio
-    scalarField thickness(pp.nPoints());
-    scalarField minThickness(pp.nPoints());
-    scalarField expansionRatio(pp.nPoints());
-    calculateLayerThickness
-    (
-        pp,
-        meshMover.adaptPatchIDs(),
-        layerParams.expansionRatio(),
-        layerParams.finalLayerRatio(),
-        layerParams.minThickness(),
-        cellLevel,
-        patchNLayers,
-        edge0Len,
-
-        thickness,
-        minThickness,
-        expansionRatio
-    );
-
-    // Calculate wall to medial axis distance for smoothing displacement
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    pointScalarField pointMedialDist
-    (
-        IOobject
-        (
-            "pointMedialDist",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        meshMover.pMesh(),
-        dimensionedScalar("pointMedialDist", dimless, 0.0)
-    );
-
-    pointVectorField dispVec
-    (
-        IOobject
-        (
-            "dispVec",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        meshMover.pMesh(),
-        dimensionedVector("dispVec", dimless, vector::zero)
-    );
-
-    pointScalarField medialRatio
-    (
-        IOobject
-        (
-            "medialRatio",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        meshMover.pMesh(),
-        dimensionedScalar("medialRatio", dimless, 0.0)
-    );
-
-    // Setup information for medial axis smoothing. Calculates medial axis
-    // and a smoothed displacement direction.
-    // - pointMedialDist : distance to medial axis
-    // - dispVec : normalised direction of nearest displacement
-    // - medialRatio : ratio of medial distance to wall distance.
-    //   (1 at wall, 0 at medial axis)
-    medialAxisSmoothingInfo
-    (
-        meshMover,
-        layerParams.nSmoothNormals(),
-        layerParams.nSmoothSurfaceNormals(),
-        layerParams.minMedianAxisAngleCos(),
-
-        dispVec,
-        medialRatio,
-        pointMedialDist
-    );
-
-
-
-    // Saved old points
-    pointField oldPoints(mesh_.points());
-
-    // Last set of topology changes. (changing mesh clears out polyTopoChange)
-    polyTopoChange savedMeshMod(mesh_.boundaryMesh().size());
-
-    boolList flaggedCells;
-    boolList flaggedFaces;
-
-    while (true)
-    {
-        // Make sure displacement is equal on both sides of coupled patches.
-        syncPatchDisplacement
-        (
-            meshMover,
-            minThickness,
-            patchDisp,
-            patchNLayers,
-            extrudeStatus
-        );
-
-        // Displacement acc. to pointnormals
-        getPatchDisplacement
-        (
-            meshMover,
-            thickness,
-            minThickness,
-            patchDisp,
-            patchNLayers,
-            extrudeStatus
-        );
-
-        // Shrink mesh by displacement value first.
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-        {
-            pointField oldPatchPos(pp.localPoints());
-
-            //// Laplacian displacement shrinking.
-            //shrinkMeshDistance
-            //(
-            //    debug,
-            //    meshMover,
-            //    -patchDisp,     // Shrink in opposite direction of addedPoints
-            //    layerParams.nSmoothDisp(),
-            //    layerParams.nSnap()
-            //);
-
-            // Medial axis based shrinking
-            shrinkMeshMedialDistance
-            (
-                meshMover,
-
-                layerParams.nSmoothThickness(),
-                layerParams.maxThicknessToMedialRatio(),
-                nAllowableErrors,
-                layerParams.nSnap(),
-                layerParams.layerTerminationCos(),
-
-                thickness,
-                minThickness,
-
-                dispVec,
-                medialRatio,
-                pointMedialDist,
-
-                extrudeStatus,
-                patchDisp,
-                patchNLayers
-            );
-
-            // Update patchDisp (since not all might have been honoured)
-            patchDisp = oldPatchPos - pp.localPoints();
-        }
-
-        // Truncate displacements that are too small (this will do internal
-        // ones, coupled ones have already been truncated by syncPatch)
-        faceSet dummySet(mesh_, "wrongPatchFaces", 0);
-        truncateDisplacement
-        (
-            meshMover,
-            minThickness,
-            dummySet,
-            patchDisp,
-            patchNLayers,
-            extrudeStatus
-        );
-
-
-        // Dump to .obj file for debugging.
-        if (debug_)
-        {
-            dumpDisplacement
-            (
-                mesh_.time().path()/"layer",
-                pp,
-                patchDisp,
-                extrudeStatus
-            );
-
-            const_cast<Time&>(mesh_.time())++;
-            Info<< "Writing shrunk mesh to " << mesh_.time().timeName() << endl;
-            mesh_.write();
-        }
-
-
-        // Mesh topo change engine
-        polyTopoChange meshMod(mesh_);
-
-        // Grow layer of cells on to patch. Handles zero sized displacement.
-        addPatchCellLayer addLayer(mesh_);
-
-        // Determine per point/per face number of layers to extrude. Also
-        // handles the slow termination of layers when going switching layers
-
-        labelList nPatchPointLayers(pp.nPoints(),-1);
-        labelList nPatchFaceLayers(pp.localFaces().size(),-1);
-        setupLayerInfoTruncation
-        (
-            meshMover,
-            patchNLayers,
-            extrudeStatus,
-            layerParams.nBufferCellsNoExtrude(),
-            nPatchPointLayers,
-            nPatchFaceLayers
-        );
-
-        // Calculate displacement for first layer for addPatchLayer
-        vectorField firstDisp(patchNLayers.size(), vector::zero);
-
-        forAll(patchNLayers, i)
-        {
-            if (patchNLayers[i] > 0)
-            {
-                if (expansionRatio[i] == 1.0)
-                {
-                    firstDisp[i] = patchDisp[i]/nPatchPointLayers[i];
-                }
-                else
-                {
-                    label nLay = nPatchPointLayers[i];
-                    scalar h =
-                        pow(expansionRatio[i], nLay - 1)
-                      * (mag(patchDisp[i])*(1.0 - expansionRatio[i]))
-                      / (1.0 - pow(expansionRatio[i], nLay));
-                    firstDisp[i] = h/mag(patchDisp[i])*patchDisp[i];
-                }
-            }
-        }
-
-        scalarField invExpansionRatio = 1.0 / expansionRatio;
-
-        // Add topo regardless of whether extrudeStatus is extruderemove.
-        // Not add layer if patchDisp is zero.
-        addLayer.setRefinement
-        (
-            invExpansionRatio,
-            pp,
-            nPatchFaceLayers,   // layers per face
-            nPatchPointLayers,  // layers per point
-            firstDisp,          // thickness of first layer
-            meshMod
-        );
-
-        if (debug_)
-        {
-            const_cast<Time&>(mesh_.time())++;
-        }
-
-        // Store mesh changes for if mesh is correct.
-        savedMeshMod = meshMod;
-
-
-        // With the stored topo changes we create a new mesh so we can
-        // undo if neccesary.
-
-        autoPtr<fvMesh> newMeshPtr;
-        autoPtr<mapPolyMesh> map = meshMod.makeMesh
-        (
-            newMeshPtr,
-            IOobject
-            (
-                //mesh.name()+"_layer",
-                mesh_.name(),
-                static_cast<polyMesh&>(mesh_).instance(),
-                mesh_.db(),
-                static_cast<polyMesh&>(mesh_).readOpt(),
-                static_cast<polyMesh&>(mesh_).writeOpt()
-            ),          // io params from original mesh but new name
-            mesh_,      // original mesh
-            true        // parallel sync
-        );
-        fvMesh& newMesh = newMeshPtr();
-
-        //?neccesary? Update fields
-        newMesh.updateMesh(map);
-
-        // Update numbering on addLayer:
-        // - cell/point labels to be newMesh.
-        // - patchFaces to remain in oldMesh order.
-        addLayer.updateMesh
-        (
-            map,
-            identity(pp.size()),
-            identity(pp.nPoints())
-        );
-
-        // Collect layer faces and cells for outside loop.
-        getLayerCellsFaces
-        (
-            newMesh,
-            addLayer,
-            flaggedCells,
-            flaggedFaces
-        );
-
-
-        if (debug_)
-        {
-            Info<< "Writing layer mesh to " << mesh_.time().timeName() << endl;
-            newMesh.write();
-            cellSet addedCellSet
-            (
-                newMesh,
-                "addedCells",
-                findIndices(flaggedCells, true)
-            );
-            Info<< "Writing "
-                << returnReduce(addedCellSet.size(), sumOp<label>())
-                << " added cells to cellSet "
-                << addedCellSet.objectPath()
-                << endl;
-            addedCellSet.write();
-
-            faceSet layerFacesSet
-            (
-                newMesh,
-                "layerFaces",
-                findIndices(flaggedCells, true)
-            );
-            Info<< "Writing "
-                << returnReduce(layerFacesSet.size(), sumOp<label>())
-                << " faces inside added layer to faceSet "
-                << layerFacesSet.objectPath()
-                << endl;
-            layerFacesSet.write();
-        }
-
-        // Unset the extrusion at the pp.
-        label nTotChanged = checkAndUnmark
-        (
-            addLayer,
-            motionDict,
-            pp,
-            newMesh,
-
-            patchDisp,
-            patchNLayers,
-            extrudeStatus
-        );
-
-        Info<< "Extruding " << countExtrusion(pp, extrudeStatus)
-            << " out of " << returnReduce(pp.size(), sumOp<label>())
-            << " faces. Removed extrusion at " << nTotChanged << " faces."
-            << endl;
-
-        if (nTotChanged == 0)
-        {
-            break;
-        }
-
-        // Reset mesh points and start again
-        meshMover.movePoints(oldPoints);
-        meshMover.correct();
-    }
-
-
-    // At this point we have a (shrunk) mesh and a set of topology changes
-    // which will make a valid mesh with layer. Apply these changes to the
-    // current mesh.
-
-    // Apply the stored topo changes to the current mesh.
-    autoPtr<mapPolyMesh> map = savedMeshMod.changeMesh(mesh_, false);
-
-    // Update fields
-    mesh_.updateMesh(map);
-
-    // Move mesh (since morphing does not do this)
-    if (map().hasMotionPoints())
-    {
-        mesh_.movePoints(map().preMotionPoints());
-    }
-
-    meshRefinerPtr_().updateMesh(map, labelList(0));
-
-
-    // Do final balancing
-    // ~~~~~~~~~~~~~~~~~~
-
-    if (Pstream::parRun())
-    {
-        // Balance. No restriction on face zones and baffles.
-        autoPtr<mapDistributePolyMesh> map = balance(false, false);
-
-        // Re-distribute flag of layer faces and cells
-        map().distributeCellData(flaggedCells);
-        map().distributeFaceData(flaggedFaces);
-    }
-
-
-    // Write mesh
-    // ~~~~~~~~~~
-
-    //writeMesh("Layer mesh");
-    cellSet addedCellSet(mesh_, "addedCells", findIndices(flaggedCells, true));
-    Info<< "Writing "
-        << returnReduce(addedCellSet.size(), sumOp<label>())
-        << " added cells to cellSet "
-        << addedCellSet.objectPath()
-        << endl;
-    addedCellSet.write();
-
-    faceSet layerFacesSet(mesh_, "layerFaces", findIndices(flaggedFaces, true));
-    Info<< "Writing "
-        << returnReduce(layerFacesSet.size(), sumOp<label>())
-        << " faces inside added layer to faceSet "
-        << layerFacesSet.objectPath()
-        << endl;
-    layerFacesSet.write();
-}
-
-
-// ************************************************************************* //
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverShrink.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverShrink.C
deleted file mode 100644
index 0c66073e84..0000000000
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverShrink.C
+++ /dev/null
@@ -1,1151 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Description
-    Shrinking mesh (part of adding cell layers)
-
-\*----------------------------------------------------------------------------*/
-
-#include "autoHexMeshDriver.H"
-#include "fvMesh.H"
-#include "Time.H"
-#include "pointFields.H"
-#include "motionSmoother.H"
-#include "pointData.H"
-#include "PointEdgeWave.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-// Calculate inverse sum of edge weights (currently always 1.0)
-void Foam::autoHexMeshDriver::sumWeights
-(
-    const PackedList<1>& isMasterEdge,
-    const labelList& meshEdges,
-    const labelList& meshPoints,
-    const edgeList& edges,
-    scalarField& invSumWeight
-) const
-{
-    invSumWeight = 0;
-
-    forAll(edges, edgeI)
-    {
-        if (isMasterEdge.get(meshEdges[edgeI]) == 1)
-        {
-            const edge& e = edges[edgeI];
-            //scalar eWeight = edgeWeights[edgeI];
-            scalar eWeight = 1.0;
-
-            invSumWeight[e[0]] += eWeight;
-            invSumWeight[e[1]] += eWeight;
-        }
-    }
-
-    syncTools::syncPointList
-    (
-        mesh_,
-        meshPoints,
-        invSumWeight,
-        plusEqOp<scalar>(),
-        scalar(0.0),        // null value
-        false               // no separation
-    );
-
-    forAll(invSumWeight, pointI)
-    {
-        scalar w = invSumWeight[pointI];
-
-        if (w > 0.0)
-        {
-            invSumWeight[pointI] = 1.0/w;
-        }
-    }
-}
-
-
-// Smooth field on moving patch
-void Foam::autoHexMeshDriver::smoothField
-(
-    const motionSmoother& meshMover,
-    const PackedList<1>& isMasterEdge,
-    const labelList& meshEdges,
-    const scalarField& fieldMin,
-    const label& nSmoothDisp,
-    scalarField& field
-) const
-{
-    const indirectPrimitivePatch& pp = meshMover.patch();
-    const edgeList& edges = pp.edges();
-    const labelList& meshPoints = pp.meshPoints();
-
-    scalarField invSumWeight(pp.nPoints());
-    sumWeights
-    (
-        isMasterEdge,
-        meshEdges,
-        meshPoints,
-        edges,
-        invSumWeight
-    );
-
-    // Get smoothly varying patch field.
-    Info<< "shrinkMeshDistance : Smoothing field ..." << endl;
-
-    for (label iter = 0; iter < nSmoothDisp; iter++)
-    {
-        scalarField average(pp.nPoints());
-        averageNeighbours
-        (
-            isMasterEdge,
-            meshEdges,
-            meshPoints,
-            pp.edges(),
-            invSumWeight,
-            field,
-            average
-        );
-
-        // Transfer to field
-        forAll(field, pointI)
-        {
-            //full smoothing neighbours + point value
-            average[pointI] = 0.5*(field[pointI]+average[pointI]);
-
-            // perform monotonic smoothing
-            if
-            (
-                average[pointI] < field[pointI]
-             && average[pointI] >= fieldMin[pointI]
-            )
-            {
-                field[pointI] = average[pointI];
-            }
-        }
-
-        // Do residual calculation every so often.
-        if ((iter % 10) == 0)
-        {
-            Info<< "    Iteration " << iter << "   residual "
-                <<  gSum(mag(field-average))
-                   /returnReduce(average.size(), sumOp<label>())
-                << endl;
-        }
-    }
-}
-
-
-// Smooth normals on moving patch.
-void Foam::autoHexMeshDriver::smoothPatchNormals
-(
-    const motionSmoother& meshMover,
-    const PackedList<1>& isMasterEdge,
-    const labelList& meshEdges,
-    const label nSmoothDisp,
-    pointField& normals
-) const
-{
-    const indirectPrimitivePatch& pp = meshMover.patch();
-    const edgeList& edges = pp.edges();
-    const labelList& meshPoints = pp.meshPoints();
-
-    // Calculate inverse sum of weights
-
-    scalarField invSumWeight(pp.nPoints());
-    sumWeights
-    (
-        isMasterEdge,
-        meshEdges,
-        meshPoints,
-        edges,
-        invSumWeight
-    );
-
-    // Get smoothly varying internal normals field.
-    Info<< "shrinkMeshDistance : Smoothing normals ..." << endl;
-
-    for (label iter = 0; iter < nSmoothDisp; iter++)
-    {
-        vectorField average(pp.nPoints());
-        averageNeighbours
-        (
-            isMasterEdge,
-            meshEdges,
-            meshPoints,
-            pp.edges(),
-            invSumWeight,
-            normals,
-            average
-        );
-
-        // Do residual calculation every so often.
-        if ((iter % 10) == 0)
-        {
-            Info<< "    Iteration " << iter << "   residual "
-                <<  gSum(mag(normals-average))
-                   /returnReduce(average.size(), sumOp<label>())
-                << endl;
-        }
-
-        // Transfer to normals vector field
-        forAll(average, pointI)
-        {
-            // full smoothing neighbours + point value
-            average[pointI] = 0.5*(normals[pointI]+average[pointI]);
-            normals[pointI] = average[pointI];
-            normals[pointI] /= mag(normals[pointI]) + VSMALL;
-        }
-    }
-}
-
-
-// Smooth normals in interior.
-void Foam::autoHexMeshDriver::smoothNormals
-(
-    const label nSmoothDisp,
-    const PackedList<1>& isMasterEdge,
-    const labelList& fixedPoints,
-    pointVectorField& normals
-) const
-{
-    // Get smoothly varying internal normals field.
-    Info<< "shrinkMeshDistance : Smoothing normals ..." << endl;
-
-    const edgeList& edges = mesh_.edges();
-
-    // Points that do not change.
-    PackedList<1> isFixedPoint(mesh_.nPoints(), 0);
-
-    // Internal points that are fixed
-    forAll(fixedPoints, i)
-    {
-        label meshPointI = fixedPoints[i];
-        isFixedPoint.set(meshPointI, 1);
-    }
-
-    // Correspondence between local edges/points and mesh edges/points
-    const labelList meshEdges(identity(mesh_.nEdges()));
-    const labelList meshPoints(identity(mesh_.nPoints()));
-
-    // Calculate inverse sum of weights
-
-    scalarField invSumWeight(mesh_.nPoints(), 0);
-    sumWeights
-    (
-        isMasterEdge,
-        meshEdges,
-        meshPoints,
-        edges,
-        invSumWeight
-    );
-
-    Info<< "shrinkMeshDistance : Smoothing normals in interior ..." << endl;
-
-    for (label iter = 0; iter < nSmoothDisp; iter++)
-    {
-        vectorField average(mesh_.nPoints());
-        averageNeighbours
-        (
-            isMasterEdge,
-            meshEdges,
-            meshPoints,
-            edges,
-            invSumWeight,
-            normals,
-            average
-        );
-
-        // Do residual calculation every so often.
-        if ((iter % 10) == 0)
-        {
-            Info<< "    Iteration " << iter << "   residual "
-                <<  gSum(mag(normals-average))
-                   /returnReduce(average.size(), sumOp<label>())
-                << endl;
-        }
-
-
-        // Transfer to normals vector field
-        forAll(average, pointI)
-        {
-            if (isFixedPoint.get(pointI) == 0)
-            {
-                //full smoothing neighbours + point value
-                average[pointI] = 0.5*(normals[pointI]+average[pointI]);
-                normals[pointI] = average[pointI];
-                normals[pointI] /= mag(normals[pointI]) + VSMALL;
-            }
-        }
-    }
-}
-
-
-// Tries and find a medial axis point. Done by comparing vectors to nearest
-// wall point for both vertices of edge.
-bool Foam::autoHexMeshDriver::isMaxEdge
-(
-    const List<pointData>& pointWallDist,
-    const label edgeI,
-    const scalar minCos
-) const
-{
-    const pointField& points = mesh_.points();
-
-    // Do not mark edges with one side on moving wall.
-
-    const edge& e = mesh_.edges()[edgeI];
-
-    vector v0(points[e[0]] - pointWallDist[e[0]].origin());
-    scalar magV0(mag(v0));
-
-    if (magV0 < SMALL)
-    {
-        return false;
-    }
-
-    vector v1(points[e[1]] - pointWallDist[e[1]].origin());
-    scalar magV1(mag(v1));
-
-    if (magV1 < SMALL)
-    {
-        return false;
-    }
-
-    v0 /= magV0;
-    v1 /= magV1;
-
-    // Test angle.
-    if ((v0 & v1) < minCos)
-    {
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
-// Stop layer growth where mesh wraps around edge with a
-// large feature angle
-void Foam::autoHexMeshDriver::handleFeatureAngleLayerTerminations
-(
-    const indirectPrimitivePatch& pp,
-    const scalar minCos,
-    List<extrudeMode>& extrudeStatus,
-    pointField& patchDisp,
-    labelList& patchNLayers,
-    label& nPointCounter
-) const
-{
-    // Mark faces that have all points extruded
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    boolList extrudedFaces(pp.size(), true);
-
-    forAll(pp.localFaces(), faceI)
-    {
-        const face& f = pp.localFaces()[faceI];
-
-        forAll(f, fp)
-        {
-            if (extrudeStatus[f[fp]] == NOEXTRUDE)
-            {
-                extrudedFaces[faceI] = false;
-                break;
-            }
-        }
-    }
-
-
-    // Detect situation where two featureedge-neighbouring faces are partly or
-    // not extruded and the edge itself is extruded. In this case unmark the
-    // edge for extrusion.
-
-    forAll(pp.edgeFaces(), edgeI)
-    {
-        const labelList& eFaces = pp.edgeFaces()[edgeI];
-
-        if (eFaces.size() == 2)
-        {
-            const edge& e = pp.edges()[edgeI];
-            label v0 = e[0];
-            label v1 = e[1];
-
-            if
-            (
-                extrudeStatus[v0] != NOEXTRUDE
-             || extrudeStatus[v1] != NOEXTRUDE
-            )
-            {
-                if (!extrudedFaces[eFaces[0]] || !extrudedFaces[eFaces[1]])
-                {
-                    const vector& n0 = pp.faceNormals()[eFaces[0]];
-                    const vector& n1 = pp.faceNormals()[eFaces[1]];
-
-                    if ((n0 & n1) < minCos)
-                    {
-                        if
-                        (
-                            unmarkExtrusion
-                            (
-                                v0,
-                                patchDisp,
-                                patchNLayers,
-                                extrudeStatus
-                            )
-                        )
-                        {
-                            nPointCounter++;
-                        }
-                        if
-                        (
-                            unmarkExtrusion
-                            (
-                                v1,
-                                patchDisp,
-                                patchNLayers,
-                                extrudeStatus
-                            )
-                        )
-                        {
-                            nPointCounter++;
-                        }
-                    }
-                }
-            }
-        }
-    }
-}
-
-
-// Find isolated islands (points, edges and faces and layer terminations)
-// in the layer mesh and stop any layer growth at these points.
-void Foam::autoHexMeshDriver::findIsolatedRegions
-(
-    const indirectPrimitivePatch& pp,
-    const PackedList<1>& isMasterEdge,
-    const labelList& meshEdges,
-    const scalar minCosLayerTermination,
-    scalarField& field,
-    List<extrudeMode>& extrudeStatus,
-    pointField& patchDisp,
-    labelList& patchNLayers
-) const
-{
-    Info<< "shrinkMeshDistance : Removing isolated regions ..." << endl;
-
-    // Keep count of number of points unextruded
-    label nPointCounter = 0;
-
-    while(true)
-    {
-        // Stop layer growth where mesh wraps around edge with a
-        // large feature angle
-        handleFeatureAngleLayerTerminations
-        (
-            pp,
-            minCosLayerTermination,
-
-            extrudeStatus,
-            patchDisp,
-            patchNLayers,
-            nPointCounter
-        );
-
-
-        // Do not extrude from point where all neighbouring
-        // faces are not grown
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-        boolList extrudedFaces(pp.size(), true);
-        forAll(pp.localFaces(), faceI)
-        {
-            const face& f = pp.localFaces()[faceI];
-            forAll(f, fp)
-            {
-                if (extrudeStatus[f[fp]] == NOEXTRUDE)
-                {
-                    extrudedFaces[faceI] = false;
-                    break;
-                }
-            }
-        }
-
-        const labelListList& pointFaces = pp.pointFaces();
-
-        boolList keptPoints(pp.nPoints(), false);
-        forAll(keptPoints, patchPointI)
-        {
-            const labelList& pFaces = pointFaces[patchPointI];
-
-            forAll(pFaces, i)
-            {
-                label faceI = pFaces[i];
-                if (extrudedFaces[faceI])
-                {
-                    keptPoints[patchPointI] = true;
-                    break;
-                }
-            }
-        }
-
-        syncTools::syncPointList
-        (
-            mesh_,
-            pp.meshPoints(),
-            keptPoints,
-            orEqOp<bool>(),
-            false,              // null value
-            false               // no separation
-        );
-
-        label nChanged = 0;
-
-        forAll(keptPoints, patchPointI)
-        {
-            if (!keptPoints[patchPointI])
-            {
-                if
-                (
-                    unmarkExtrusion
-                    (
-                        patchPointI,
-                        patchDisp,
-                        patchNLayers,
-                        extrudeStatus
-                    )
-                )
-                {
-                   nPointCounter++;
-                   nChanged++;
-                }
-           }
-       }
-
-
-       if (returnReduce(nChanged, sumOp<label>()) == 0)
-       {
-           break;
-       }
-   }
-
-    const edgeList& edges = pp.edges();
-
-
-    // Count number of mesh edges using a point
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    labelList isolatedPoint(pp.nPoints(),0);
-
-    forAll(edges, edgeI)
-    {
-        if (isMasterEdge.get(meshEdges[edgeI]) == 1)
-        {
-            const edge& e = edges[edgeI];
-
-            label v0 = e[0];
-            label v1 = e[1];
-
-            if (extrudeStatus[v1] != NOEXTRUDE)
-            {
-                isolatedPoint[v0] += 1;
-            }
-            if (extrudeStatus[v0] != NOEXTRUDE)
-            {
-                isolatedPoint[v1] += 1;
-            }
-        }
-    }
-
-    syncTools::syncPointList
-    (
-        mesh_,
-        pp.meshPoints(),
-        isolatedPoint,
-        plusEqOp<label>(),
-        0,       // null value
-        false    // no separation
-    );
-
-    // stop layer growth on isolated faces
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    forAll(pp, faceI)
-    {
-        const face& f = pp.localFaces()[faceI];
-        bool failed = false;
-        forAll(f, fp)
-        {
-            if (isolatedPoint[f[fp]] > 2)
-            {
-                failed = true;
-                break;
-            }
-        }
-        bool allPointsExtruded = true;
-        if (!failed)
-        {
-            forAll(f, fp)
-            {
-                if (extrudeStatus[f[fp]] == NOEXTRUDE)
-                {
-                    allPointsExtruded = false;
-                    break;
-                }
-            }
-
-            if (allPointsExtruded)
-            {
-                forAll(f, fp)
-                {
-                    if
-                    (
-                        unmarkExtrusion
-                        (
-                            f[fp],
-                            patchDisp,
-                            patchNLayers,
-                            extrudeStatus
-                        )
-                    )
-                    {
-                        field[f[fp]] = 0.0;
-                        nPointCounter++;
-                    }
-                }
-            }
-        }
-    }
-
-    reduce(nPointCounter, sumOp<label>());
-    Info<< "Number isolated points extrusion stopped : "<< nPointCounter<< endl;
-}
-
-
-// Calculates medial axis fields:
-// dispVec     : normal of nearest wall. Where this normal changes direction
-//               defines the medial axis
-// medialDist  : distance to medial axis
-// medialRatio : ratio of medial distance to wall distance.
-//               (1 at wall, 0 at medial axis)
-void Foam::autoHexMeshDriver::medialAxisSmoothingInfo
-(
-    const motionSmoother& meshMover,
-    const label nSmoothNormals,
-    const label nSmoothSurfaceNormals,
-    const scalar minMedianAxisAngleCos,
-
-    pointVectorField& dispVec,
-    pointScalarField& medialRatio,
-    pointScalarField& medialDist
-) const
-{
-
-    Info<< "medialAxisSmoothingInfo :"
-        << " Calculate distance to Medial Axis ..." << endl;
-
-    const pointField& points = mesh_.points();
-    const pointMesh& pMesh = meshMover.pMesh();
-
-    const indirectPrimitivePatch& pp = meshMover.patch();
-    const vectorField& faceNormals = pp.faceNormals();
-    const labelList& meshPoints = pp.meshPoints();
-
-    // Predetermine mesh edges
-    // ~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Precalulate master edge (only relevant for shared edges)
-    PackedList<1> isMasterEdge(syncTools::getMasterEdges(mesh_));
-    // Precalculate meshEdge per pp edge
-    labelList meshEdges(pp.nEdges());
-
-    forAll(meshEdges, patchEdgeI)
-    {
-        const edge& e = pp.edges()[patchEdgeI];
-
-        label v0 = pp.meshPoints()[e[0]];
-        label v1 = pp.meshPoints()[e[1]];
-        meshEdges[patchEdgeI] = meshTools::findEdge
-        (
-            mesh_.edges(),
-            mesh_.pointEdges()[v0],
-            v0,
-            v1
-        );
-    }
-
-
-    // Determine pointNormal
-    // ~~~~~~~~~~~~~~~~~~~~~
-
-    pointField pointNormals(pp.nPoints(), vector::zero);
-    {
-        labelList nPointFaces(pp.nPoints(), 0);
-
-        forAll(faceNormals, faceI)
-        {
-            const face& f = pp.localFaces()[faceI];
-
-            forAll(f, fp)
-            {
-                pointNormals[f[fp]] += faceNormals[faceI];
-                nPointFaces[f[fp]] ++;
-            }
-        }
-
-        syncTools::syncPointList
-        (
-            mesh_,
-            meshPoints,
-            pointNormals,
-            plusEqOp<vector>(),
-            vector::zero,       // null value
-            false               // no separation
-        );
-
-        syncTools::syncPointList
-        (
-            mesh_,
-            meshPoints,
-            nPointFaces,
-            plusEqOp<label>(),
-            0,                  // null value
-            false               // no separation
-        );
-
-        forAll(pointNormals, i)
-        {
-            pointNormals[i] /= nPointFaces[i];
-        }
-    }
-
-    // Smooth patch normal vectors
-    smoothPatchNormals
-    (
-        meshMover,
-        isMasterEdge,
-        meshEdges,
-        nSmoothSurfaceNormals,
-        pointNormals
-    );
-
-
-    // Calculate distance to pp points
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Distance to wall
-    List<pointData> pointWallDist(mesh_.nPoints());
-
-
-    // 1. Calculate distance to points where displacement is specified.
-    {
-        // Seed data.
-        List<pointData> wallInfo(meshPoints.size());
-
-        forAll(meshPoints, patchPointI)
-        {
-            label pointI = meshPoints[patchPointI];
-            wallInfo[patchPointI] = pointData
-            (
-                points[pointI],
-                0.0,
-                pointI,                       // passive scalar
-                pointNormals[patchPointI]     // surface normals
-            );
-        }
-
-        // Do all calculations
-        List<pointData> edgeWallDist(mesh_.nEdges());
-        PointEdgeWave<pointData> wallDistCalc
-        (
-            pMesh,
-            meshPoints,
-            wallInfo,
-            pointWallDist,
-            edgeWallDist,
-            mesh_.nPoints()  // max iterations
-        );
-    }
-
-    // 2. Find points with max distance and transport information back to
-    //    wall.
-    {
-        List<pointData> pointMedialDist(mesh_.nPoints());
-        List<pointData> edgeMedialDist(mesh_.nEdges());
-
-        // Seed point data.
-        DynamicList<pointData> maxInfo(meshPoints.size());
-        DynamicList<label> maxPoints(meshPoints.size());
-
-        // 1. Medial axis points
-
-        const edgeList& edges = mesh_.edges();
-
-        forAll(edges, edgeI)
-        {
-            if (isMaxEdge(pointWallDist, edgeI, minMedianAxisAngleCos))
-            {
-                // Both end points of edge have very different nearest wall
-                // point. Mark both points as medial axis points.
-                const edge& e = edges[edgeI];
-
-                forAll(e, ep)
-                {
-                    label pointI = e[ep];
-
-                    if (!pointMedialDist[pointI].valid())
-                    {
-                        maxPoints.append(pointI);
-                        maxInfo.append
-                        (
-                            pointData
-                            (
-                                points[pointI],
-                                0.0,
-                                pointI,         // passive data
-                                vector::zero    // passive data
-                            )
-                        );
-                        pointMedialDist[pointI] = maxInfo[maxInfo.size()-1];
-                    }
-                }
-            }
-        }
-
-
-        // 2. Seed non-adapt patches
-        const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-
-        labelHashSet adaptPatches(meshMover.adaptPatchIDs());
-
-        forAll(patches, patchI)
-        {
-            const polyPatch& pp = patches[patchI];
-
-            if
-            (
-                !pp.coupled()
-             && !isA<emptyPolyPatch>(pp)
-             && !adaptPatches.found(patchI)
-            )
-            {
-                Info<< "Inserting points on patch " << pp.name() << endl;
-
-                const labelList& meshPoints = pp.meshPoints();
-
-                forAll(meshPoints, i)
-                {
-                    label pointI = meshPoints[i];
-
-                    if (!pointMedialDist[pointI].valid())
-                    {
-                        maxPoints.append(pointI);
-                        maxInfo.append
-                        (
-                            pointData
-                            (
-                                points[pointI],
-                                0.0,
-                                pointI,         // passive data
-                                vector::zero    // passive data
-                            )
-                        );
-                        pointMedialDist[pointI] = maxInfo[maxInfo.size()-1];
-                    }
-                }
-            }
-        }
-
-        maxInfo.shrink();
-        maxPoints.shrink();
-
-        // Do all calculations
-        PointEdgeWave<pointData> medialDistCalc
-        (
-            pMesh,
-            maxPoints,
-            maxInfo,
-
-            pointMedialDist,
-            edgeMedialDist,
-            mesh_.nPoints()  // max iterations
-        );
-
-        // Extract medial axis distance as pointScalarField
-        forAll(pointMedialDist, pointI)
-        {
-            medialDist[pointI] = Foam::sqrt(pointMedialDist[pointI].distSqr());
-        }
-    }
-
-    // Extract transported surface normals as pointVectorField
-    forAll(dispVec, i)
-    {
-        dispVec[i] = pointWallDist[i].v();
-    }
-
-    // Smooth normal vectors. Do not change normals on pp.meshPoints
-    smoothNormals(nSmoothNormals, isMasterEdge, meshPoints, dispVec);
-
-    // Calculate ratio point medial distance to point wall distance
-    forAll(medialRatio, pointI)
-    {
-        scalar wDist2 = pointWallDist[pointI].distSqr();
-        scalar mDist = medialDist[pointI];
-
-        if (wDist2 < sqr(SMALL) && mDist < SMALL)
-        {
-            medialRatio[pointI] = 0.0;
-        }
-        else
-        {
-            medialRatio[pointI] = mDist / (Foam::sqrt(wDist2) + mDist);
-        }
-    }
-
-    if (debug_)
-    {
-        Info<< "medialAxisSmoothingInfo :"
-            << " Writing:" << nl
-            << "    " << dispVec.name()
-            << " : normalised direction of nearest displacement" << nl
-            << "    " << medialDist.name()
-            << " : distance to medial axis" << nl
-            << "    " << medialRatio.name()
-            << " : ratio of medial distance to wall distance" << nl
-            << endl;
-        dispVec.write();
-        medialDist.write();
-        medialRatio.write();
-    }
-}
-
-
-void Foam::autoHexMeshDriver::shrinkMeshMedialDistance
-(
-    motionSmoother& meshMover,
-    const label nSmoothThickness,
-    const scalar maxThicknessToMedialRatio,
-    const label nAllowableErrors,
-    const label nSnap,
-    const scalar minCosLayerTermination,
-
-    const scalarField& layerThickness,
-    const scalarField& minThickness,
-
-    const pointVectorField& dispVec,
-    const pointScalarField& medialRatio,
-    const pointScalarField& medialDist,
-
-    List<extrudeMode>& extrudeStatus,
-    pointField& patchDisp,
-    labelList& patchNLayers
-) const
-{
-    Info<< "shrinkMeshMedialDistance : Smoothing using Medial Axis ..." << endl;
-
-    const polyMesh& mesh = meshMover.mesh();
-    const pointMesh& pMesh = meshMover.pMesh();
-
-    const indirectPrimitivePatch& pp = meshMover.patch();
-    const labelList& meshPoints = pp.meshPoints();
-
-    // Precalulate master edge (only relevant for shared edges)
-    PackedList<1> isMasterEdge(syncTools::getMasterEdges(mesh_));
-    // Precalculate meshEdge per pp edge
-    labelList meshEdges(pp.nEdges());
-
-    forAll(meshEdges, patchEdgeI)
-    {
-        const edge& e = pp.edges()[patchEdgeI];
-
-        label v0 = pp.meshPoints()[e[0]];
-        label v1 = pp.meshPoints()[e[1]];
-        meshEdges[patchEdgeI] = meshTools::findEdge
-        (
-            mesh_.edges(),
-            mesh_.pointEdges()[v0],
-            v0,
-            v1
-        );
-    }
-
-
-    scalarField thickness(layerThickness.size());
-
-    thickness = mag(patchDisp);
-
-    forAll(thickness, patchPointI)
-    {
-        if (extrudeStatus[patchPointI] == NOEXTRUDE)
-        {
-            thickness[patchPointI] = 0.0;
-        }
-    }
-
-    label numThicknessRatioExclude = 0;
-
-    // reduce thickness where thickness/medial axis distance large
-    forAll(meshPoints, patchPointI)
-    {
-        if (extrudeStatus[patchPointI] != NOEXTRUDE)
-        {
-            label pointI = meshPoints[patchPointI];
-
-            scalar mDist = medialDist[pointI];
-
-            scalar thicknessRatio = thickness[patchPointI]/(mDist+VSMALL);
-
-            if (thicknessRatio > maxThicknessToMedialRatio)
-            {
-                // Truncate thickness.
-                thickness[patchPointI] =
-                    0.5*(minThickness[patchPointI]+thickness[patchPointI]);
-
-                patchDisp[patchPointI] =
-                    thickness[patchPointI]
-                  * patchDisp[patchPointI]
-                  / (mag(patchDisp[patchPointI]) + VSMALL);
-                numThicknessRatioExclude++;
-            }
-        }
-    }
-
-    reduce(numThicknessRatioExclude, sumOp<label>());
-    Info<< "shrinkMeshMedialDistance : Reduce layer thickness at "
-        << numThicknessRatioExclude
-        << " nodes where thickness to medial axis distance is large " << endl;
-
-    // find points where layer growth isolated to a lone point, edge or face
-
-    findIsolatedRegions
-    (
-        pp,
-        isMasterEdge,
-        meshEdges,
-        minCosLayerTermination,
-        thickness,
-        extrudeStatus,
-        patchDisp,
-        patchNLayers
-    );
-
-    // smooth layer thickness on moving patch
-    smoothField
-    (
-        meshMover,
-        isMasterEdge,
-        meshEdges,
-        minThickness,
-        nSmoothThickness,
-        thickness
-    );
-
-    List<pointData> pointWallDist(mesh.nPoints());
-
-    const pointField& points = mesh.points();
-    // 1. Calculate distance to points where displacement is specified.
-    // This wave is used to transport layer thickness
-    {
-        // Distance to wall and medial axis on edges.
-        List<pointData> edgeWallDist(mesh.nEdges());
-        labelList wallPoints(meshPoints.size());
-
-        // Seed data.
-        List<pointData> wallInfo(meshPoints.size());
-
-        forAll(meshPoints, patchPointI)
-        {
-            label pointI = meshPoints[patchPointI];
-            wallPoints[patchPointI] = pointI;
-            wallInfo[patchPointI] = pointData
-            (
-                points[pointI],
-                0.0,
-                thickness[patchPointI],       // transport layer thickness
-                vector::zero                  // passive vector
-            );
-        }
-
-        // Do all calculations
-        PointEdgeWave<pointData> wallDistCalc
-        (
-            pMesh,
-            wallPoints,
-            wallInfo,
-            pointWallDist,
-            edgeWallDist,
-            mesh.nPoints()  // max iterations
-        );
-    }
-
-    // Calculate scaled displacement vector
-    pointVectorField& displacement = meshMover.displacement();
-
-    forAll(displacement, pointI)
-    {
-        // 1) displacement on nearest wall point, scaled by medialRatio
-        //    (wall distance / medial distance)
-        // 2) pointWallDist[pointI].s() is layer thickness transported
-        //    from closest wall point.
-        // 3) shrink in opposite direction of addedPoints
-        displacement[pointI] =
-            -medialRatio[pointI]
-            *pointWallDist[pointI].s()
-            *dispVec[pointI];
-    }
-
-    // Current faces to check. Gets modified in meshMover.scaleMesh
-    labelList checkFaces(identity(mesh.nFaces()));
-
-    Info<< "shrinkMeshMedialDistance : Moving mesh ..." << endl;
-    scalar oldErrorReduction = -1;
-
-    for (label iter = 0; iter < 2*nSnap ; iter++)
-    {
-        Info<< "Iteration " << iter << endl;
-        if (iter == nSnap)
-        {
-            Info<< "Displacement scaling for error reduction set to 0." << endl;
-            oldErrorReduction = meshMover.setErrorReduction(0.0);
-        }
-
-        if (meshMover.scaleMesh(checkFaces, true, nAllowableErrors))
-        {
-            Info<< "shrinkMeshMedialDistance : Successfully moved mesh" << endl;
-            break;
-        }
-    }
-
-    if (oldErrorReduction >= 0)
-    {
-        meshMover.setErrorReduction(oldErrorReduction);
-    }
-
-    Info<< "shrinkMeshMedialDistance : Finished moving mesh ..." << endl;
-}
-
-
-// ************************************************************************* //
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverSnap.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverSnap.C
deleted file mode 100644
index 8e074cd05e..0000000000
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverSnap.C
+++ /dev/null
@@ -1,1245 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Description
-    All to do with snapping to the surface
-
-\*----------------------------------------------------------------------------*/
-
-#include "autoHexMeshDriver.H"
-#include "syncTools.H"
-#include "fvMesh.H"
-#include "Time.H"
-#include "OFstream.H"
-#include "mapPolyMesh.H"
-#include "motionSmoother.H"
-#include "pointEdgePoint.H"
-#include "PointEdgeWave.H"
-#include "mergePoints.H"
-#include "snapParameters.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-void Foam::autoHexMeshDriver::getZonedSurfaces
-(
-    labelList& zonedSurfaces,
-    labelList& unzonedSurfaces
-) const
-{    // Surfaces with zone information
-    const wordList& faceZoneNames = surfaces().faceZoneNames();
-
-    zonedSurfaces.setSize(faceZoneNames.size());
-    label zonedI = 0;
-    unzonedSurfaces.setSize(faceZoneNames.size());
-    label unzonedI = 0;
-
-    forAll(faceZoneNames, surfI)
-    {
-        if (faceZoneNames[surfI].size() > 0)
-        {
-            zonedSurfaces[zonedI++] = surfI;
-        }
-        else
-        {
-            unzonedSurfaces[unzonedI++] = surfI;
-        }
-    }
-    zonedSurfaces.setSize(zonedI);
-    unzonedSurfaces.setSize(unzonedI);
-}
-
-
-// Get faces to repatch. Returns map from face to patch.
-Foam::Map<Foam::label> Foam::autoHexMeshDriver::getZoneBafflePatches
-(
-    const bool allowBoundary
-) const
-{
-    Map<label> bafflePatch(mesh_.nFaces()/1000);
-
-    const wordList& faceZoneNames = surfaces().faceZoneNames();
-    const faceZoneMesh& fZones = mesh_.faceZones();
-
-    forAll(faceZoneNames, surfI)
-    {
-        if (faceZoneNames[surfI].size() > 0)
-        {
-            // Get zone
-            label zoneI = fZones.findZoneID(faceZoneNames[surfI]);
-
-            const faceZone& fZone = fZones[zoneI];
-
-            //// Get patch allocated for zone
-            //label patchI = surfaceToCyclicPatch_[surfI];
-            // Get patch of (first region) of surface
-            label patchI = globalToPatch_[surfaces().globalRegion(surfI, 0)];
-
-            Info<< "For surface "
-                << surfaces().names()[surfI]
-                << " found faceZone " << fZone.name()
-                << " and patch " << mesh_.boundaryMesh()[patchI].name()
-                << endl;
-
-
-            forAll(fZone, i)
-            {
-                label faceI = fZone[i];
-
-                if (allowBoundary || mesh_.isInternalFace(faceI))
-                {
-                    if (!bafflePatch.insert(faceI, patchI))
-                    {
-                        label oldPatchI = bafflePatch[faceI];
-
-                        if (oldPatchI != patchI)
-                        {
-                            FatalErrorIn("getZoneBafflePatches(const bool)")
-                                << "Face " << faceI
-                                << " fc:" << mesh_.faceCentres()[faceI]
-                                << " is in faceZone "
-                                << mesh_.boundaryMesh()[oldPatchI].name()
-                                << " and in faceZone "
-                                << mesh_.boundaryMesh()[patchI].name()
-                                << abort(FatalError);
-                        }
-                    }
-                }
-            }
-        }
-    }
-    return bafflePatch;
-}
-
-
-// Calculate geometrically collocated points, Requires PackedList to be
-// sizes and initalised!
-Foam::label Foam::autoHexMeshDriver::getCollocatedPoints
-(
-    const scalar tol,
-    const pointField& points,
-    PackedList<1>& isCollocatedPoint
-)
-{
-    labelList pointMap;
-    pointField newPoints;
-    bool hasMerged = mergePoints
-    (
-        points,                         // points
-        tol,                            // mergeTol
-        false,                          // verbose
-        pointMap,
-        newPoints
-    );
-
-    if (!returnReduce(hasMerged, orOp<bool>()))
-    {
-        return 0;
-    }
-
-    // Determine which newPoints 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(newPoints.size(), -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>());
-}
-
-
-// Calculate displacement as average of patch points.
-Foam::pointField Foam::autoHexMeshDriver::smoothPatchDisplacement
-(
-    const motionSmoother& meshMover
-) const
-{
-    const indirectPrimitivePatch& pp = meshMover.patch();
-
-    // Calculate geometrically non-manifold points on the patch to be moved.
-    PackedList<1> nonManifoldPoint(pp.nPoints());
-    label nNonManifoldPoints = getCollocatedPoints
-    (
-        SMALL,
-        pp.localPoints(),
-        nonManifoldPoint
-    );
-    Info<< "Found " << nNonManifoldPoints << " non-mainfold 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 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)
-        {
-            avgBoundary[patchPointI] += pp[pFaces[pfI]].centre(points);
-        }
-        nBoundary[patchPointI] = pFaces.size();
-    }
-
-    syncTools::syncPointList
-    (
-        mesh,
-        pp.meshPoints(),
-        avgBoundary,
-        plusEqOp<point>(),  // combine op
-        vector::zero,       // null value
-        false               // no separation
-    );
-    syncTools::syncPointList
-    (
-        mesh,
-        pp.meshPoints(),
-        nBoundary,
-        plusEqOp<label>(),  // combine op
-        0,                  // null value
-        false               // no separation
-    );
-
-    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 (Pstream::parRun() && isA<processorPolyPatch>(patches[patchI]))
-            {
-                const processorPolyPatch& pp =
-                    refCast<const processorPolyPatch>(patches[patchI]);
-
-                if (pp.myProcNo() < pp.neighbProcNo())
-                {
-                    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]]++;
-                        }
-                    }
-                }
-            }
-            else if (isA<cyclicPolyPatch>(patches[patchI]))
-            {
-                const cyclicPolyPatch& pp =
-                    refCast<const cyclicPolyPatch>(patches[patchI]);
-
-                const vectorField::subField faceCentres = pp.faceCentres();
-
-                for (label i = 0; i < pp.size()/2; 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
-            false               // no separation
-        );
-        syncTools::syncPointList
-        (
-            mesh,
-            globalNum,
-            plusEqOp<label>(),  // combine op
-            0,                  // null value
-            false               // no separation
-        );
-
-        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.faceNeighbour(), faceI)
-    {
-        label own = mesh.faceOwner()[faceI];
-        const face& f = mesh.faces()[faceI];
-
-        forAll(f, fp)
-        {
-            anyCell[f[fp]] = own;
-        }
-    }
-    for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
-    {
-        label own = mesh.faceOwner()[faceI];
-
-        const face& f = mesh.faces()[faceI];
-
-        forAll(f, fp)
-        {
-            anyCell[f[fp]] = own;
-        }
-    }
-
-
-    // Displacement to calculate.
-    pointField patchDisp(meshPoints.size(), vector::zero);
-
-    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) == 0u)
-        {
-            // 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 patchDisp;
-}
-
-
-Foam::tmp<Foam::scalarField> Foam::autoHexMeshDriver::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
-    (
-        pMesh,
-        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());
-    }
-
-
-    //{
-    //    // For debugging: dump to file
-    //    pointScalarField pointDist
-    //    (
-    //        IOobject
-    //        (
-    //            "pointDist",
-    //            mesh.DB().timeName(),
-    //            mesh.DB(),
-    //            IOobject::NO_READ,
-    //            IOobject::AUTO_WRITE
-    //        ),
-    //        pMesh,
-    //        dimensionedScalar("pointDist", dimless, 0.0)
-    //    );
-    //
-    //    forAll(allEdgeInfo, edgeI)
-    //    {
-    //        scalar d = Foam::sqrt(allEdgeInfo[edgeI].distSqr());
-    //
-    //        const edge& e = mesh.edges()[edgeI];
-    //
-    //        pointDist[e[0]] += d;
-    //        pointDist[e[1]] += d;
-    //    }
-    //    forAll(pointDist, pointI)
-    //    {
-    //        pointDist[pointI] /= mesh.pointEdges()[pointI].size();
-    //    }
-    //    Info<< "Writing patch distance to " << pointDist.name()
-    //        << " at time " << mesh.DB().timeName() << endl;
-    //
-    //    pointDist.write();
-    //}
-
-    return tedgeDist;
-}
-
-
-void Foam::autoHexMeshDriver::dumpMove
-(
-    const fileName& fName,
-    const pointField& meshPts,
-    const pointField& surfPts
-)
-{
-    // Dump direction of growth into file
-    Pout<< nl << "Dumping move direction to " << fName << nl
-        << "View this Lightwave-OBJ file with e.g. javaview" << nl
-        << 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::autoHexMeshDriver::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;
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-Foam::autoPtr<Foam::mapPolyMesh> Foam::autoHexMeshDriver::createZoneBaffles
-(
-    List<labelPair>& baffles
-)
-{
-    labelList zonedSurfaces;
-    labelList unzonedSurfaces;
-    getZonedSurfaces(zonedSurfaces, unzonedSurfaces);
-
-    autoPtr<mapPolyMesh> map;
-
-    // No need to sync; all processors will have all same zonedSurfaces.
-    if (zonedSurfaces.size() > 0)
-    {
-        // Split internal faces on interface surfaces
-        Info<< "Converting zoned faces into baffles ..." << endl;
-
-        // Get faces (internal only) to be baffled. Map from face to patch
-        // label.
-        Map<label> faceToPatch(getZoneBafflePatches(false));
-
-        label nZoneFaces = returnReduce(faceToPatch.size(), sumOp<label>());
-        if (nZoneFaces > 0)
-        {
-            // Convert into labelLists
-            labelList ownPatch(mesh_.nFaces(), -1);
-            forAllConstIter(Map<label>, faceToPatch, iter)
-            {
-                ownPatch[iter.key()] = iter();
-            }
-
-            // Create baffles. both sides same patch.
-            map = meshRefinerPtr_().createBaffles(ownPatch, ownPatch);
-
-            // Get pairs of faces created.
-            // Just loop over faceMap and store baffle if we encounter a slave
-            // face.
-
-            baffles.setSize(faceToPatch.size());
-            label baffleI = 0;
-
-            const labelList& faceMap = map().faceMap();
-            const labelList& reverseFaceMap = map().reverseFaceMap();
-
-            forAll(faceMap, faceI)
-            {
-                label oldFaceI = faceMap[faceI];
-
-                // Does face originate from face-to-patch
-                Map<label>::const_iterator iter = faceToPatch.find(oldFaceI);
-
-                if (iter != faceToPatch.end())
-                {
-                    label masterFaceI = reverseFaceMap[oldFaceI];
-                    if (faceI != masterFaceI)
-                    {
-                        baffles[baffleI++] = labelPair(masterFaceI, faceI);
-                    }
-                }
-            }
-
-            if (baffleI != faceToPatch.size())
-            {
-                FatalErrorIn("autoHexMeshDriver::createZoneBaffles(..)")
-                    << "Had " << faceToPatch.size() << " patches to create "
-                    << " but encountered " << baffleI
-                    << " slave faces originating from patcheable faces."
-                    << abort(FatalError);
-            }
-
-            if (debug_)
-            {
-                const_cast<Time&>(mesh_.time())++;
-                Pout<< "Writing baffled mesh to time "
-                    << mesh_.time().timeName() << endl;
-                mesh_.write();
-            }
-        }
-        Info<< "Created " << nZoneFaces << " baffles in = "
-            << mesh_.time().cpuTimeIncrement() << " s\n" << nl << endl;
-    }
-    return map;
-}
-
-
-Foam::autoPtr<Foam::mapPolyMesh> Foam::autoHexMeshDriver::mergeZoneBaffles
-(
-    const List<labelPair>& baffles
-)
-{
-    labelList zonedSurfaces;
-    labelList unzonedSurfaces;
-    getZonedSurfaces(zonedSurfaces, unzonedSurfaces);
-
-    autoPtr<mapPolyMesh> map;
-
-    // No need to sync; all processors will have all same zonedSurfaces.
-    label nBaffles = returnReduce(baffles.size(), sumOp<label>());
-    if (zonedSurfaces.size() > 0 && nBaffles > 0)
-    {
-        // Merge any baffles
-        Info<< "Converting " << nBaffles << " baffles back into zoned faces ..."
-            << endl;
-
-        map = meshRefinerPtr_().mergeBaffles(baffles);
-
-        Info<< "Converted baffles in = "
-            << mesh_.time().cpuTimeIncrement() << " s\n" << nl << endl;
-    }
-
-    return map;
-}
-
-
-Foam::scalarField Foam::autoHexMeshDriver::calcSnapDistance
-(
-    const snapParameters& snapParams,
-    const indirectPrimitivePatch& pp
-) const
-{
-    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
-        false               // no separation
-    );
-
-    return snapParams.snapTol()*maxEdgeLen;
-}
-
-
-// Invert globalToPatch_ to get the patches related to surfaces.
-Foam::labelList Foam::autoHexMeshDriver::getSurfacePatches() const
-{
-    // Set of patches originating from surface
-    labelHashSet surfacePatchSet(globalToPatch_.size());
-
-    forAll(globalToPatch_, i)
-    {
-        if (globalToPatch_[i] != -1)
-        {
-            surfacePatchSet.insert(globalToPatch_[i]);
-        }
-    }
-
-    DynamicList<label> surfacePatches(surfacePatchSet.size());
-
-    for (label patchI = 0; patchI < mesh_.boundaryMesh().size(); patchI++)
-    {
-        if (surfacePatchSet.found(patchI))
-        {
-            surfacePatches.append(patchI);
-        }
-    }
-    return surfacePatches.shrink();
-}
-
-
-void Foam::autoHexMeshDriver::preSmoothPatch
-(
-    const snapParameters& snapParams,
-    const label nInitErrors,
-    const List<labelPair>& baffles,
-    motionSmoother& meshMover
-) const
-{
-    labelList checkFaces;
-
-    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;
-        }
-
-        pointField patchDisp(smoothPatchDisplacement(meshMover));
-
-        // The current mesh is the starting mesh to smooth from.
-        meshMover.setDisplacement(patchDisp);
-        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_)
-    {
-        const_cast<Time&>(mesh_.time())++;
-        Pout<< "Writing patch smoothed mesh to time " << mesh_.time().timeName()
-            << endl;
-        mesh_.write();
-    }
-
-    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::autoHexMeshDriver::getZoneSurfacePoints
-(
-    const indirectPrimitivePatch& pp,
-    const word& zoneName
-) const
-{
-    label zoneI = mesh_.faceZones().findZoneID(zoneName);
-
-    if (zoneI == -1)
-    {
-        FatalErrorIn
-        (
-            "autoHexMeshDriver::getZoneSurfacePoints"
-            "(const indirectPrimitivePatch&, const word&)"
-        )   << "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::vectorField Foam::autoHexMeshDriver::calcNearestSurface
-(
-    const scalarField& snapDist,
-    motionSmoother& meshMover
-) const
-{
-    Info<< "Calculating patchDisplacement as distance to nearest surface"
-        << " point ..." << endl;
-
-    const indirectPrimitivePatch& pp = meshMover.patch();
-    const pointField& localPoints = pp.localPoints();
-
-    // Divide surfaces into zoned and unzoned
-    labelList zonedSurfaces;
-    labelList unzonedSurfaces;
-    getZonedSurfaces(zonedSurfaces, unzonedSurfaces);
-
-    // Displacement per patch point
-    vectorField patchDisp(localPoints.size(), vector::zero);
-
-
-    // 1. All points to non-interface surfaces
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    {
-        List<pointIndexHit> hitInfo;
-        labelList hitSurface;
-        surfaces().findNearest
-        (
-            unzonedSurfaces,
-            localPoints,
-            sqr(4*snapDist),        // sqr of attract distance
-            hitSurface,
-            hitInfo
-        );
-
-        forAll(hitInfo, pointI)
-        {
-            if (hitInfo[pointI].hit())
-            {
-                patchDisp[pointI] =
-                    hitInfo[pointI].hitPoint()
-                  - localPoints[pointI];
-            }
-            //else
-            //{
-            //   WarningIn("autoHexMeshDriver::calcNearestSurface(..)")
-            //        << "For point:" << pointI
-            //        << " coordinate:" << localPoints[pointI]
-            //        << " did not find any surface within:"
-            //        << 4*snapDist[pointI]
-            //        << " meter." << endl;
-            //}
-        }
-    }
-
-
-
-    // 2. All points on zones to their respective surface
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Surfaces with zone information
-    const wordList& faceZoneNames = surfaces().faceZoneNames();
-
-    forAll(zonedSurfaces, i)
-    {
-        label zoneSurfI = zonedSurfaces[i];
-
-        const labelList surfacesToTest(1, zoneSurfI);
-
-        // Get indices of points both on faceZone and on pp.
-        labelList zonePointIndices
-        (
-            getZoneSurfacePoints
-            (
-                pp,
-                faceZoneNames[zoneSurfI]
-            )
-        );
-
-        pointField zonePoints(zonePointIndices.size());
-        forAll(zonePointIndices, i)
-        {
-            zonePoints[i] = localPoints[zonePointIndices[i]];
-        }
-
-        // Find nearest for points both on faceZone and pp.
-        List<pointIndexHit> hitInfo;
-        labelList hitSurface;
-        surfaces().findNearest
-        (
-            labelList(1, zoneSurfI),
-            zonePoints,
-            sqr(4*snapDist),
-            hitSurface,
-            hitInfo
-        );
-
-        forAll(hitInfo, pointI)
-        {
-            if (hitInfo[pointI].hit())
-            {
-                patchDisp[pointI] =
-                    hitInfo[pointI].hitPoint()
-                  - localPoints[pointI];
-            }
-            else
-            {
-                WarningIn("autoHexMeshDriver::calcNearestSurface(..)")
-                    << "For point:" << pointI
-                    << " coordinate:" << localPoints[pointI]
-                    << " did not find any surface within:"
-                    << 4*snapDist[pointI]
-                    << " meter." << endl;
-            }
-        }
-    }
-
-
-    {
-        scalarField magDisp(mag(patchDisp));
-
-        Info<< "Wanted displacement : average:"
-            << gSum(magDisp)/returnReduce(patchDisp.size(), sumOp<label>())
-            << " min:" << gMin(magDisp)
-            << " max:" << gMax(magDisp) << endl;
-    }
-
-    Info<< "Calculated surface displacement in = "
-        << mesh_.time().cpuTimeIncrement() << " s\n" << nl << endl;
-
-
-    // Limit amount of movement.
-    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,
-        minMagEqOp(),                   // combine op
-        vector(GREAT, GREAT, GREAT),    // null value
-        false                           // no separation
-    );
-
-
-    // Check for displacement being outwards.
-    outwardsDisplacement(pp, patchDisp);
-
-    // Set initial distribution of displacement field (on patches) from
-    // patchDisp and make displacement consistent with b.c. on displacement
-    // pointVectorField.
-    meshMover.setDisplacement(patchDisp);
-
-    if (debug_)
-    {
-        dumpMove
-        (
-            mesh_.time().path()/"patchDisplacement.obj",
-            pp.localPoints(),
-            pp.localPoints() + patchDisp
-        );
-    }
-
-    return patchDisp;
-}
-
-
-void Foam::autoHexMeshDriver::smoothDisplacement
-(
-    const snapParameters& snapParams,
-    motionSmoother& meshMover
-) const
-{
-    const pointMesh& pMesh = meshMover.pMesh();
-    const indirectPrimitivePatch& pp = meshMover.patch();
-
-    Info<< "Smoothing displacement ..." << endl;
-
-    // Set edge diffusivity as inverse of distance to patch
-    scalarField edgeGamma(1.0/(edgePatchDist(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, false, disp);
-    }
-    Info<< "Displacement smoothed in = "
-        << mesh_.time().cpuTimeIncrement() << " s\n" << nl << endl;
-
-    if (debug_)
-    {
-        const_cast<Time&>(mesh_.time())++;
-        Pout<< "Writing smoothed mesh to time " << mesh_.time().timeName()
-            << endl;
-        mesh_.write();
-
-        Pout<< "Writing displacement field ..." << endl;
-        disp.write();
-        tmp<pointScalarField> magDisp(mag(disp));
-        magDisp().write();
-
-        Pout<< "Writing actual patch displacement ..." << endl;
-        vectorField actualPatchDisp
-        (
-            IndirectList<point>(disp, pp.meshPoints())()
-        );
-        dumpMove
-        (
-            mesh_.time().path()/"actualPatchDisplacement.obj",
-            pp.localPoints(),
-            pp.localPoints() + actualPatchDisp
-        );
-    }
-}
-
-
-void Foam::autoHexMeshDriver::scaleMesh
-(
-    const snapParameters& snapParams,
-    const label nInitErrors,
-    const List<labelPair>& baffles,
-    motionSmoother& meshMover
-)
-{
-    // Relax displacement until correct mesh
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    labelList checkFaces(identity(mesh_.nFaces()));
-
-    scalar oldErrorReduction = -1;
-
-    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);
-        }
-
-        if (meshMover.scaleMesh(checkFaces, baffles, true, nInitErrors))
-        {
-            Info<< "Successfully moved mesh" << endl;
-
-            break;
-        }
-        if (debug_)
-        {
-            const_cast<Time&>(mesh_.time())++;
-            Pout<< "Writing scaled mesh to time " << mesh_.time().timeName()
-                << endl;
-            mesh_.write();
-
-            Pout<< "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;
-}
-
-
-// ************************************************************************* //
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverTemplates.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverTemplates.C
deleted file mode 100644
index 31807f0dee..0000000000
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverTemplates.C
+++ /dev/null
@@ -1,75 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "autoHexMeshDriver.H"
-#include "syncTools.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-template<class Type>
-void Foam::autoHexMeshDriver::averageNeighbours
-(
-    const PackedList<1>& isMasterEdge,
-    const labelList& meshEdges,
-    const labelList& meshPoints,
-    const edgeList& edges,
-    const scalarField& invSumWeight,
-    const Field<Type>& data,
-    Field<Type>& average
-) const
-{
-    average = pTraits<Type>::zero;
-
-    forAll(edges, edgeI)
-    {
-        if (isMasterEdge.get(meshEdges[edgeI]) == 1)
-        {
-            const edge& e = edges[edgeI];
-            //scalar eWeight = edgeWeights[edgeI];
-            scalar eWeight =  1.0;
-            label v0 = e[0];
-            label v1 = e[1];
-
-            average[v0] += eWeight*data[v1];
-            average[v1] += eWeight*data[v0];
-        }
-    }
-
-    syncTools::syncPointList
-    (
-        mesh_,
-        meshPoints,
-        average,
-        plusEqOp<Type>(),
-        pTraits<Type>::zero,    // null value
-        false                   // no separation
-    );
-
-    average *= invSumWeight;
-}
-
-
-// ************************************************************************* //
diff --git a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.H b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.H
index d764c3fa6f..2db82c8eed 100644
--- a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.H
+++ b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.H
@@ -361,7 +361,7 @@ public:
 
     // Member Functions
 
-        //- Helper funntion: count cells per processor in wanted distribution
+        //- Helper function: count cells per processor in wanted distribution
         static labelList countCells(const labelList&);
 
         //- Send cells to neighbours according to distribution