ENH: snappyHexMesh: various improvements. See below or the default snappyHexMeshDict.

Refinement:
-----------
// Optionally avoid patch merging - keeps hexahedral cells
// (to be used with automatic refinement/unrefinement)
//mergePatchFaces off;

// Optional multiple locationsInMesh with corresponding optional cellZone
// (automatically generates faceZones inbetween)
locationsInMesh
(
    ((-0.09 -0.039 -0.049)  bottomAir)  // cellZone bottomAir
    ((-0.09 0.009 -0.049)   topAir)     // cellZone topAir
);

// Optional faceType and patchType specification for these faceZones
faceZoneControls
{
    bottomAir_to_topAir
    {
        faceType baffle;
    }
}

/ Optional checking of 'bleeding' of mesh through a specifying a locations
// outside the mesh
locationsOutsideMesh ((0 0 0)(12.3 101.17 3.98));

// Improved refinement: refine all cells with all (or all but one) sides refined

// Improved refinement: refine all cells with opposing faces with different
// refinement level. These cells can happen on multiply curved surfaces.
// Default on, can be switched off with
//interfaceRefine false;

Snapping
--------
// Optional smoothing of points at refinement interfaces. This will reduce
// the non-orthogonality at refinement interfaces.
//nSmoothInternal $nSmoothPatch;

Layering
--------

// Layers can be added to patches or to any side of a faceZone.
// (Any faceZone internally gets represented as two patches)

// The angle to merge patch faces can be set independently of the
// featureAngle. This is especially useful for large feature angles
// Default is the same as the featureAngle.
//mergePatchFacesAngle 45;

// Optional mesh shrinking type 'displacementMotionSolver'. It uses any
// displacementMotionSolver, e.g. displacementSBRStress
// (default is the medial-axis algorithm, 'displacementMedialAxis')
//meshShrinker displacementMotionSolver;

applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -493,21 +493,30 @@ int main(int argc, char *argv[])
         //    new patchID to neighbour processor)
         //  - number of new patches (nPatches)
 
-        labelList sidePatchID;
+        labelList edgePatchID;
+        labelList edgeZoneID;
+        boolList edgeFlip;
+        labelList inflateFaceID;
         label nPatches;
         Map<label> nbrProcToPatch;
         Map<label> patchToNbrProc;
-        addPatchCellLayer::calcSidePatch
+        addPatchCellLayer::calcExtrudeInfo
         (
+            true,   // for internal edges get zone info from any face
+
             mesh,
             globalFaces,
             edgeGlobalFaces,
             extrudePatch,
 
-            sidePatchID,
+            edgePatchID,
             nPatches,
             nbrProcToPatch,
-            patchToNbrProc
+            patchToNbrProc,
+
+            edgeZoneID,
+            edgeFlip,
+            inflateFaceID
         );
 
 
@@ -659,7 +668,12 @@ int main(int argc, char *argv[])
 
             ratio,              // expansion ratio
             extrudePatch,       // patch faces to extrude
-            sidePatchID,        // if boundary edge: patch to use
+
+            edgePatchID,        // if boundary edge: patch for extruded face
+            edgeZoneID,         // optional zone for extruded face
+            edgeFlip,
+            inflateFaceID,      // mesh face that zone/patch info is from
+
             exposedPatchID,     // if new mesh: patches for exposed faces
             nFaceLayers,
             nPointLayers,

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -1121,6 +1121,38 @@ int main(int argc, char *argv[])
     );
 
 
+    // Refinement parameters
+    const refinementParameters refineParams(refineDict);
+
+    // Snap parameters
+    const snapParameters snapParams(snapDict);
+
+
+
+    // Add all the cellZones and faceZones
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // 1. cellZones relating to surface (faceZones added later)
+
+    const labelList namedSurfaces
+    (
+        surfaceZonesInfo::getNamedSurfaces(surfaces.surfZones())
+    );
+
+    labelList surfaceToCellZone = surfaceZonesInfo::addCellZonesToMesh
+    (
+        surfaces.surfZones(),
+        namedSurfaces,
+        mesh
+    );
+
+
+    // 2. cellZones relating to locations
+
+    refineParams.addCellZonesToMesh(mesh);
+
+
+
     // Add all the surface regions as patches
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -1128,6 +1160,8 @@ int main(int argc, char *argv[])
     //  (faceZone surfaces)
     labelList globalToMasterPatch;
     labelList globalToSlavePatch;
+
+
     {
         Info<< nl
             << "Adding patches for surface regions" << nl
@@ -1148,6 +1182,7 @@ int main(int argc, char *argv[])
 
         const labelList& surfaceGeometry = surfaces.surfaces();
         const PtrList<dictionary>& surfacePatchInfo = surfaces.patchInfo();
+        const polyBoundaryMesh& pbm = mesh.boundaryMesh();
 
         forAll(surfaceGeometry, surfI)
         {
@@ -1157,7 +1192,9 @@ int main(int argc, char *argv[])
 
             Info<< surfaces.names()[surfI] << ':' << nl << nl;
 
-            if (surfaces.surfZones()[surfI].faceZoneName().empty())
+            const word& fzName = surfaces.surfZones()[surfI].faceZoneName();
+
+            if (fzName.empty())
             {
                 // 'Normal' surface
                 forAll(regNames, i)
@@ -1188,7 +1225,7 @@ int main(int argc, char *argv[])
 
                     Info<< setf(ios_base::left)
                         << setw(6) << patchI
-                        << setw(20) << mesh.boundaryMesh()[patchI].type()
+                        << setw(20) << pbm[patchI].type()
                         << setw(30) << regNames[i] << nl;
 
                     globalToMasterPatch[globalRegionI] = patchI;
@@ -1228,7 +1265,7 @@ int main(int argc, char *argv[])
 
                         Info<< setf(ios_base::left)
                             << setw(6) << patchI
-                            << setw(20) << mesh.boundaryMesh()[patchI].type()
+                            << setw(20) << pbm[patchI].type()
                             << setw(30) << regNames[i] << nl;
 
                         globalToMasterPatch[globalRegionI] = patchI;
@@ -1260,12 +1297,27 @@ int main(int argc, char *argv[])
 
                         Info<< setf(ios_base::left)
                             << setw(6) << patchI
-                            << setw(20) << mesh.boundaryMesh()[patchI].type()
+                            << setw(20) << pbm[patchI].type()
                             << setw(30) << slaveName << nl;
 
                         globalToSlavePatch[globalRegionI] = patchI;
                     }
                 }
+
+                // For now: have single faceZone per surface. Use first
+                // region in surface for patch for zoneing
+                if (regNames.size())
+                {
+                    label globalRegionI = surfaces.globalRegion(surfI, 0);
+
+                    meshRefiner.addFaceZone
+                    (
+                        fzName,
+                        pbm[globalToMasterPatch[globalRegionI]].name(),
+                        pbm[globalToSlavePatch[globalRegionI]].name(),
+                        surfaces.surfZones()[surfI].faceType()
+                    );
+                }
             }
 
             Info<< nl;
@@ -1275,10 +1327,85 @@ int main(int argc, char *argv[])
     }
 
 
+
+    // Add all information for all the remaining faceZones
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    HashTable<Pair<word> > faceZoneToPatches;
+    forAll(mesh.faceZones(), zoneI)
+    {
+        const word& fzName = mesh.faceZones()[zoneI].name();
+
+        label mpI, spI;
+        surfaceZonesInfo::faceZoneType fzType;
+        bool hasInfo = meshRefiner.getFaceZoneInfo(fzName, mpI, spI, fzType);
+
+        if (!hasInfo)
+        {
+            // faceZone does not originate from a surface but presumably
+            // from a cellZone pair instead
+            string::size_type i = fzName.find("_to_");
+            if (i != string::npos)
+            {
+                word cz0 = fzName.substr(0, i);
+                word cz1 = fzName.substr(i+4, fzName.size()-i+4);
+                word slaveName(cz1 + "_to_" + cz0);
+                faceZoneToPatches.insert(fzName, Pair<word>(fzName, slaveName));
+            }
+            else
+            {
+                // Add as fzName + fzName_slave
+                const word slaveName = fzName + "_slave";
+                faceZoneToPatches.insert(fzName, Pair<word>(fzName, slaveName));
+            }
+        }
+    }
+
+    if (faceZoneToPatches.size())
+    {
+        autoRefineDriver::addFaceZones
+        (
+            meshRefiner,
+            refineParams,
+            faceZoneToPatches
+        );
+    }
+
+
+
+    // Re-do intersections on meshed boundaries since they use an extrapolated
+    // other side
+    {
+        const labelList adaptPatchIDs(meshRefiner.meshedPatches());
+
+        const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+
+        label nFaces = 0;
+        forAll(adaptPatchIDs, i)
+        {
+            nFaces += pbm[adaptPatchIDs[i]].size();
+        }
+
+        labelList faceLabels(nFaces);
+        nFaces = 0;
+        forAll(adaptPatchIDs, i)
+        {
+            const polyPatch& pp = pbm[adaptPatchIDs[i]];
+            forAll(pp, i)
+            {
+                faceLabels[nFaces++] = pp.start()+i;
+            }
+        }
+        meshRefiner.updateIntersections(faceLabels);
+    }
+
+
+
     // Parallel
     // ~~~~~~~~
 
-    // Decomposition
+    // Construct decomposition engine. Note: cannot use decompositionModel
+    // MeshObject since we're clearing out the mesh inside the mesh generation.
     autoPtr<decompositionMethod> decomposerPtr
     (
         decompositionMethod::New
@@ -1312,14 +1439,17 @@ int main(int argc, char *argv[])
     const Switch wantSnap(meshDict.lookup("snap"));
     const Switch wantLayers(meshDict.lookup("addLayers"));
 
-    // Refinement parameters
-    const refinementParameters refineParams(refineDict);
+    const Switch mergePatchFaces
+    (
+        meshDict.lookupOrDefault("mergePatchFaces", true)
+    );
 
-    // Snap parameters
-    const snapParameters snapParams(snapDict);
-
-    // Layer addition parameters
-    const layerParameters layerParams(layerDict, mesh.boundaryMesh());
+    if (!mergePatchFaces)
+    {
+        Info<< "Not merging patch-faces of cell to preserve"
+            << " (split)hex cell shape."
+            << nl << endl;
+    }
 
 
     if (wantRefine)
@@ -1348,6 +1478,7 @@ int main(int argc, char *argv[])
             refineParams,
             snapParams,
             refineParams.handleSnapProblems(),
+            mergePatchFaces,        // merge co-planar faces
             motionDict
         );
 
@@ -1387,6 +1518,7 @@ int main(int argc, char *argv[])
         (
             snapDict,
             motionDict,
+            mergePatchFaces,
             curvature,
             planarAngle,
             snapParams
@@ -1408,6 +1540,9 @@ int main(int argc, char *argv[])
     {
         cpuTime timer;
 
+        // Layer addition parameters
+        const layerParameters layerParams(layerDict, mesh.boundaryMesh());
+
         autoLayerDriver layerDriver
         (
             meshRefiner,
@@ -1433,6 +1568,7 @@ int main(int argc, char *argv[])
             layerDict,
             motionDict,
             layerParams,
+            mergePatchFaces,
             preBalance,
             decomposer,
             distributor

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict

@@ -71,6 +71,11 @@ geometry
     }
 };
 
+
+// Optional: avoid patch-face merging. Allows mesh to be used for
+//           refinement/unrefinement
+//mergePatchFaces    off; // default on
+
 // Settings for the castellatedMesh generation.
 castellatedMeshControls
 {
@@ -177,7 +182,7 @@ castellatedMeshControls
             //  how to select the cells that are in the cellZone
             //  (inside / outside / specified insidePoint)
             //  The orientation of the faceZone is
-            //  - if on cellZone(s) : point out of (maximum) cellZone
+            //  - if on cellZone(s) : point out of (minimum) cellZone
             //  - if freestanding   : oriented according to surface
 
             //faceZone sphere;
@@ -249,16 +254,70 @@ castellatedMeshControls
 
     // After refinement patches get added for all refinementSurfaces and
     // all cells intersecting the surfaces get put into these patches. The
-    // section reachable from the locationInMesh is kept.
+    // section reachable from the location(s)InMesh is kept.
     // NOTE: This point should never be on a face, always inside a cell, even
     // after refinement.
-    locationInMesh (5 0.28 0.43);
+    //
+    // There are two different ways of specifying the regions to keep:
+    // 1. a single locationInMesh. All the 'zoned' surfaces are marked as such
+    //    in the refinementSurfaces with the faceZone and cellZone keywords.
+    //
+    // or
+    //
+    // 2. multiple locationsInMesh, with per location the name of the cellZone.
+    //    This uses walking to determine zones and automatically creates
+    //    faceZones on the outside of cellZones.
 
-    // Whether any faceZones (as specified in the refinementSurfaces)
-    // are only on the boundary of corresponding cellZones or also allow
-    // free-standing zone faces. Not used if there are no faceZones.
-    allowFreeStandingZoneFaces true;
 
+        // Ad 1. Specify a single location and how to treat faces inbetween
+        //       cellZones
+        locationInMesh (5 0.28 0.43);
+
+        // Whether any faceZones (as specified in the refinementSurfaces)
+        // are only on the boundary of corresponding cellZones or also allow
+        // free-standing zone faces. Not used if there are no faceZones.
+        allowFreeStandingZoneFaces true;
+
+
+
+        // 2. Specify multiple locations with optional cellZones for the
+        //    regions. faceZones are automatically constructed from the
+        //    names of the cellZones: <cellZoneA> _to_ <cellZoneB>
+        //    where the cellZoneA is the lowest numbered cellZone (non-cellZone
+        //    cells are in a special region called "none" which is always
+        //    last).
+
+        locationsInMesh
+        (
+            ((-0.09 -0.039 -0.049)  bottomAir)  // cellZone 0
+            ((-0.09 0.009 -0.049)   topAir)     // cellZone 1
+            ((-0.09 0.001 -0.049)   leftSolid)  // cellZone 2
+            ((0.02  0.001 -0.049)   rightSolid) // cellZone 3
+            ((-0.001 -0.039 0.0015) heater)     // cellZone 4
+        );
+
+        // Per synthesised faceZone name the faceType and type of baffles to
+        // create
+        faceZoneControls
+        {
+            bottomAir_to_heater
+            {
+                // Optional specification of patch type (default is wall). No
+                // constraint types (cyclic, symmetry) etc. are allowed.
+                patchInfo
+                {
+                    type patch;
+                    inGroups (patchPatches);
+                }
+                faceType baffle;
+            }
+        }
+
+
+
+    // Optional locations that should not be reachable from
+    // location(s)InMesh
+    locationsOutsideMesh ((100 100 100));
 
     // Optional: do not remove cells likely to give snapping problems
     // handleSnapProblems false;
@@ -266,6 +325,10 @@ castellatedMeshControls
     // Optional: switch off topological test for cells to-be-squashed
     //           and use geometric test instead
     //useTopologicalSnapDetection false;
+
+    // Optional: do not refine surface cells with opposite faces of
+    //           differing refinement levels
+    //interfaceRefine false;
 }
 
 // Settings for the snapping.
@@ -275,6 +338,11 @@ snapControls
     // to surface
     nSmoothPatch 3;
 
+    // Optional: number of smoothing iterations for internal points on
+    // refinement interfaces. This will reduce non-orthogonality on
+    // refinement interfaces.
+    //nSmoothInternal $nSmoothPatch;
+
     // Maximum relative distance for points to be attracted by surface.
     // True distance is this factor times local maximum edge length.
     // Note: changed(corrected) w.r.t 17x! (17x used 2* tolerance)
@@ -287,6 +355,11 @@ snapControls
     // before upon reaching a correct mesh.
     nRelaxIter 5;
 
+    // (wip) disable snapping to opposite near surfaces (revert to 22x
+    //  behaviour)
+    // detectNearSurfacesSnap false;
+
+
     // Feature snapping
 
         // Number of feature edge snapping iterations.
@@ -305,8 +378,28 @@ snapControls
         multiRegionFeatureSnap false;
 
 
-    // wip: disable snapping to opposite near surfaces (revert to 22x behaviour)
-    // detectNearSurfacesSnap false;
+        //- When to run face splitting (never at first iteration, always
+        //  at last iteration). Is interval. Default -1 (disabled)
+        //nFaceSplitInterval 5;
+
+
+        // (wip) Optional for explicit feature snapping:
+        //- Detect baffle edges. Default is true.
+        //detectBaffles false;
+        //- Erase attraction close to feature point. Default is false.
+        //releasePoints true;
+        //- Walk along feature edges, adding missing ones. Default is true.
+        //stringFeatures false;
+        //- If diagonal attraction also attract other face points. Default is
+        //  false
+        //avoidDiagonal true;
+        //- When splitting what concave faces to leave intact. Default is 45
+        //  degrees.
+        //concaveAngle 30;
+        //- When splitting the minimum area ratio of faces. If face split
+        //  causes ratio of area less than this do not split. Default is 0.3
+        //minAreaRatio 0.3;
+
 }
 
 // Settings for the layer addition.
@@ -350,10 +443,10 @@ addLayersControls
     // cannot be above minThickness do not add layer.
     // If relativeSizes this is relative to undistorted size of cell
     // outside layer..
-    minThickness 0.25;
+    minThickness 0.1;
 
 
-    // Per final patch (so not geometry!) the layer information
+    // Per final patch or faceZone (so not geometry!) the layer information
     // Note: This behaviour changed after 21x. Any non-mentioned patches
     //       now slide unless:
     //          - nSurfaceLayers is explicitly mentioned to be 0.
@@ -397,6 +490,10 @@ addLayersControls
         // are perpendicular
         featureAngle 130;
 
+        // When to merge patch faces. Default is featureAngle. Useful when
+        // featureAngle is large.
+        //mergePatchFacesAngle 45;
+
         // Stop layer growth on highly warped cells
         maxFaceThicknessRatio 0.5;
 
@@ -433,8 +530,10 @@ addLayersControls
         // default is 0.
         //nSmoothDisplacement 90;
 
-        // (wip)Optional: do not extrude a point if none of the surrounding points is
-        // not extruded. Default is false.
+        // (wip)Optional: do not extrude any point where
+        //   (false) : all surrounding faces are not fully extruded
+        //   (true)  : all surrounding points are not extruded
+        // Default is false.
         //detectExtrusionIsland true;
 
 
@@ -449,7 +548,8 @@ addLayersControls
         // before upon reaching a correct mesh.
         nRelaxIter 5;
 
-        // Create buffer region for new layer terminations
+        // Create buffer region for new layer terminations, i.e. gradually
+        // step down number of layers. Set to <0 to terminate layer in one go.
         nBufferCellsNoExtrude 0;
 
         // Overall max number of layer addition iterations. The mesher will

src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -34,6 +34,8 @@ License
 #include "polyModifyFace.H"
 #include "polyAddCell.H"
 #include "globalIndex.H"
+#include "PatchTools.H"
+#include "dummyTransform.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -221,13 +223,15 @@ Foam::labelPair Foam::addPatchCellLayer::getEdgeString
 }
 
 
-// Adds a side face i.e. extrudes a patch edge.
 Foam::label Foam::addPatchCellLayer::addSideFace
 (
     const indirectPrimitivePatch& pp,
     const labelListList& addedCells,    // per pp face the new extruded cell
     const face& newFace,
     const label newPatchID,
+    const label zoneI,
+    const bool edgeFlip,
+    const label inflateFaceI,
 
     const label ownFaceI,               // pp face that provides owner
     const label nbrFaceI,
@@ -238,68 +242,14 @@ Foam::label Foam::addPatchCellLayer::addSideFace
     polyTopoChange& meshMod
 ) const
 {
-    // Face or edge to 'inflate' from
-    label inflateEdgeI = -1;
-    label inflateFaceI = -1;
-
-    // Check mesh faces using edge
-    if (addToMesh_)
-    {
-        forAll(meshFaces, i)
-        {
-            if (mesh_.isInternalFace(meshFaces[i]))
-            {
-                // meshEdge uses internal faces so ok to inflate from it
-                inflateEdgeI = meshEdgeI;
-                break;
-            }
-        }
-    }
-
-    // Zone info comes from any side patch face. Otherwise -1 since we
-    // don't know what to put it in - inherit from the extruded faces?
-    label zoneI = -1;   //mesh_.faceZones().whichZone(meshFaceI);
-    bool flip = false;
-
     label addedFaceI = -1;
 
+
     // Is patch edge external edge of indirectPrimitivePatch?
     if (nbrFaceI == -1)
     {
         // External edge so external face.
 
-        const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-
-        // Loop over all faces connected to edge to inflate and
-        // see if we can find a face that is otherPatchID
-
-        // Get my mesh face and its zone.
-        label meshFaceI = pp.addressing()[ownFaceI];
-
-        forAll(meshFaces, k)
-        {
-            label faceI = meshFaces[k];
-
-            if
-            (
-                (faceI != meshFaceI)
-             && (patches.whichPatch(faceI) == newPatchID)
-            )
-            {
-                // Found the patch face. Use it to inflate from
-                inflateEdgeI = -1;
-                inflateFaceI = faceI;
-
-                zoneI = mesh_.faceZones().whichZone(faceI);
-                if (zoneI != -1)
-                {
-                    label index = mesh_.faceZones()[zoneI].whichFace(faceI);
-                    flip = mesh_.faceZones()[zoneI].flipMap()[index];
-                }
-                break;
-            }
-        }
-
         // Determine if different number of layer on owner and neighbour side
         // (relevant only for coupled faces). See section for internal edge
         // below.
@@ -333,16 +283,16 @@ Foam::label Foam::addPatchCellLayer::addSideFace
         (
             polyAddFace
             (
-                newFace,                    // face
-                addedCells[ownFaceI][layerOwn],   // owner
-                -1,                         // neighbour
-                -1,                         // master point
-                inflateEdgeI,               // master edge
-                inflateFaceI,               // master face
-                false,                      // flux flip
-                newPatchID,                 // patch for face
-                zoneI,                      // zone for face
-                flip                        // face zone flip
+                newFace,                        // face
+                addedCells[ownFaceI][layerOwn], // owner
+                -1,                             // neighbour
+                -1,                             // master point
+                -1,                             // master edge
+                inflateFaceI,                   // master face
+                false,                          // flux flip
+                newPatchID,                     // patch for face
+                zoneI,                          // zone for face
+                edgeFlip                        // face zone flip
             )
         );
     }
@@ -395,20 +345,37 @@ Foam::label Foam::addPatchCellLayer::addSideFace
             layerOwn = layerI;
         }
 
+
+        // Check mesh internal faces using edge to initialise
+        label inflateEdgeI = -1;
+        if (addToMesh_)
+        {
+            forAll(meshFaces, i)
+            {
+                if (mesh_.isInternalFace(meshFaces[i]))
+                {
+                    // meshEdge uses internal faces so ok to inflate from it
+                    inflateEdgeI = meshEdgeI;
+                    break;
+                }
+            }
+        }
+
+
         addedFaceI = meshMod.setAction
         (
             polyAddFace
             (
-                newFace,                    // face
-                addedCells[ownFaceI][layerOwn],   // owner
-                addedCells[nbrFaceI][layerNbr],   // neighbour
-                -1,                         // master point
-                inflateEdgeI,               // master edge
-                -1,                         // master face
-                false,                      // flux flip
-                -1,                         // patch for face
-                zoneI,                      // zone for face
-                flip                        // face zone flip
+                newFace,                        // face
+                addedCells[ownFaceI][layerOwn], // owner
+                addedCells[nbrFaceI][layerNbr], // neighbour
+                -1,                             // master point
+                inflateEdgeI,                   // master edge
+                -1,                             // master face
+                false,                          // flux flip
+                -1,                             // patch for face
+                zoneI,                          // zone for face
+                edgeFlip                        // face zone flip
             )
         );
 
@@ -467,6 +434,134 @@ void Foam::addPatchCellLayer::setFaceProps
 }
 
 
+void Foam::addPatchCellLayer::setFaceProps
+(
+    const polyMesh& mesh,
+    const indirectPrimitivePatch& pp,
+    const label ppEdgeI,
+    const label faceI,
+
+    label& patchI,
+    label& zoneI,
+    bool& zoneFlip,
+    label& inflateFaceI
+)
+{
+    setFaceProps
+    (
+        mesh,
+        faceI,
+
+        patchI,
+        zoneI,
+        zoneFlip
+    );
+
+    if (patchI != -1 || zoneI != -1)
+    {
+        inflateFaceI = faceI;
+    }
+
+    if (zoneI != -1)
+    {
+        // Correct flip for patch edge ordering
+        const edge& ppEdge = pp.edges()[ppEdgeI];
+        const edge patchEdge
+        (
+            pp.meshPoints()[ppEdge[0]],
+            pp.meshPoints()[ppEdge[1]]
+        );
+
+        const face& f = mesh.faces()[faceI];
+        bool found = false;
+        forAll(f, fp)
+        {
+            const edge e(f[fp], f.nextLabel(fp));
+            int stat = edge::compare(e, patchEdge);
+            if (stat == 1)
+            {
+                found = true;
+                break;
+            }
+            else if (stat == -1)
+            {
+                found = true;
+                zoneFlip = !zoneFlip;
+                break;
+            }
+        }
+
+        if (!found)
+        {
+            FatalErrorIn("addPatchCellLayer::setFaceProps(..)")
+                << "Problem: cannot find patch edge " << ppEdgeI
+                << " with mesh vertices " << patchEdge
+                << " at " << patchEdge.line(mesh.points())
+                << " is not in face " << faceI << " with mesh vertices "
+                << f
+                << exit(FatalError);
+        }
+    }
+}
+
+
+void Foam::addPatchCellLayer::findZoneFace
+(
+    const bool useInternalFaces,
+    const bool useBoundaryFaces,
+
+    const polyMesh& mesh,
+    const indirectPrimitivePatch& pp,
+    const label ppEdgeI,
+    const UIndirectList<label>& excludeFaces,
+    const labelList& meshFaces,
+
+    label& inflateFaceI,
+    label& patchI,
+    label& zoneI,
+    bool& zoneFlip
+)
+{
+    inflateFaceI = -1;
+    patchI = -1;
+    zoneI = -1;
+    zoneFlip = false;
+
+    forAll(meshFaces, k)
+    {
+        label faceI = meshFaces[k];
+
+        if
+        (
+            (findIndex(excludeFaces, faceI) == -1)
+         && (
+                (mesh.isInternalFace(faceI) && useInternalFaces)
+             || (!mesh.isInternalFace(faceI) && useBoundaryFaces)
+            )
+        )
+        {
+            setFaceProps
+            (
+                mesh,
+                pp,
+                ppEdgeI,
+                faceI,
+
+                patchI,
+                zoneI,
+                zoneFlip,
+                inflateFaceI
+            );
+
+            if (zoneI != -1 || patchI != -1)
+            {
+                break;
+            }
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from mesh
@@ -564,38 +659,61 @@ Foam::labelListList Foam::addPatchCellLayer::globalEdgeFaces
 }
 
 
-void Foam::addPatchCellLayer::calcSidePatch
+void Foam::addPatchCellLayer::calcExtrudeInfo
 (
+    const bool zoneFromAnyFace,
+
     const polyMesh& mesh,
     const globalIndex& globalFaces,
     const labelListList& globalEdgeFaces,
     const indirectPrimitivePatch& pp,
 
-    labelList& sidePatchID,
+    labelList& edgePatchID,
     label& nPatches,
     Map<label>& nbrProcToPatch,
-    Map<label>& patchToNbrProc
+    Map<label>& patchToNbrProc,
+    labelList& edgeZoneID,
+    boolList& edgeFlip,
+    labelList& inflateFaceID
 )
 {
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
+    const globalMeshData& gd = mesh.globalData();
 
     // Precalculate mesh edges for pp.edges.
     const labelList meshEdges(pp.meshEdges(mesh.edges(), mesh.pointEdges()));
 
-    sidePatchID.setSize(pp.nEdges());
-    sidePatchID = -1;
+    edgePatchID.setSize(pp.nEdges());
+    edgePatchID = -1;
+    edgeZoneID.setSize(pp.nEdges());
+    edgeZoneID = -1;
+    edgeFlip.setSize(pp.nEdges());
+    edgeFlip = false;
+
+
+    // Determine properties for faces extruded from edges
+    // - edge inbetween two different processors:
+    //      - extrude as patch face on correct processor
+    // - edge at end of patch (so edgeFaces.size() == 1):
+    //      - use mesh face that is a boundary face
+    // - edge internal to patch (so edgeFaces.size() == 2):
+
+
 
     // These also get determined but not (yet) exported:
     // - whether face is created from other face or edge
-    // - what zone&orientation face should have
 
-    labelList inflateEdgeI(pp.nEdges(), -1);
-    labelList inflateFaceI(pp.nEdges(), -1);
-    labelList sideZoneID(pp.nEdges(), -1);
-    boolList sideFlip(pp.nEdges(), false);
+    inflateFaceID.setSize(pp.nEdges(), -1);
 
     nPatches = patches.size();
 
+
+
+    // Pass1:
+    // For all edges inbetween two processors: see if matches to existing
+    // processor patch or create interprocessor-patch if necessary.
+    // Sets edgePatchID[edgeI] but none of the other quantities
+
     forAll(globalEdgeFaces, edgeI)
     {
         const labelList& eGlobalFaces = globalEdgeFaces[edgeI];
@@ -606,8 +724,7 @@ void Foam::addPatchCellLayer::calcSidePatch
         )
         {
             // Locally but not globally a boundary edge. Hence a coupled
-            // edge. Find the patch to use if on different
-            // processors.
+            // edge. Find the patch to use if on different processors.
 
             label f0 = eGlobalFaces[0];
             label f1 = eGlobalFaces[1];
@@ -625,18 +742,25 @@ void Foam::addPatchCellLayer::calcSidePatch
 
             if (otherProcI != -1)
             {
-                sidePatchID[edgeI] = findProcPatch(mesh, otherProcI);
-                if (sidePatchID[edgeI] == -1)
+                if (findIndex(gd[Pstream::myProcNo()], otherProcI) != -1)
+                {
+                    // There is already a processorPolyPatch to otherProcI.
+                    // Use it. Note that we can only index procPatchMap
+                    // if the processor actually is a neighbour processor
+                    // so that is why we first check.
+                    edgePatchID[edgeI] = gd.procPatchMap()[otherProcI];
+                }
+                else
                 {
                     // Cannot find a patch to processor. See if already
                     // marked for addition
                     if (nbrProcToPatch.found(otherProcI))
                     {
-                        sidePatchID[edgeI] = nbrProcToPatch[otherProcI];
+                        edgePatchID[edgeI] = nbrProcToPatch[otherProcI];
                     }
                     else
                     {
-                        sidePatchID[edgeI] = nPatches;
+                        edgePatchID[edgeI] = nPatches;
                         nbrProcToPatch.insert(otherProcI, nPatches);
                         patchToNbrProc.insert(nPatches, otherProcI);
                         nPatches++;
@@ -647,9 +771,8 @@ void Foam::addPatchCellLayer::calcSidePatch
     }
 
 
-
-    // Determine face properties for all other boundary edges
-    // ------------------------------------------------------
+    // Pass2: determine face properties for all other edges
+    // ----------------------------------------------------
 
     const labelListList& edgeFaces = pp.edgeFaces();
 
@@ -657,13 +780,10 @@ void Foam::addPatchCellLayer::calcSidePatch
 
     forAll(edgeFaces, edgeI)
     {
-        if (edgeFaces[edgeI].size() == 1 && sidePatchID[edgeI] == -1)
+        if (edgePatchID[edgeI] == -1)
         {
-            // Proper, uncoupled patch edge.
+            UIndirectList<label> ppFaces(pp.addressing(), edgeFaces[edgeI]);
 
-            label myFaceI = pp.addressing()[edgeFaces[edgeI][0]];
-
-            // Pick up any boundary face on this edge and use its properties
             label meshEdgeI = meshEdges[edgeI];
             const labelList& meshFaces = mesh.edgeFaces
             (
@@ -671,43 +791,71 @@ void Foam::addPatchCellLayer::calcSidePatch
                 dynMeshEdgeFaces
             );
 
-            forAll(meshFaces, k)
+            if (edgeFaces[edgeI].size() == 2)
             {
-                label faceI = meshFaces[k];
+                // Internal edge. Look at any face (internal or boundary) to
+                // determine extrusion properties. First one that has zone
+                // info wins
 
-                if (faceI != myFaceI && !mesh.isInternalFace(faceI))
-                {
-                    setFaceProps
-                    (
-                        mesh,
-                        faceI,
+                label dummyPatchI = -1;
+                findZoneFace
+                (
+                    true,               // useInternalFaces,
+                    zoneFromAnyFace,    // useBoundaryFaces,
 
-                        sidePatchID[edgeI],
-                        sideZoneID[edgeI],
-                        sideFlip[edgeI]
-                    );
-                    inflateFaceI[edgeI] = faceI;
-                    inflateEdgeI[edgeI] = -1;
+                    mesh,
+                    pp,
+                    edgeI,
 
-                    break;
-                }
+
+                    ppFaces,            // excludeFaces,
+                    meshFaces,          // meshFaces,
+
+                    inflateFaceID[edgeI],
+                    dummyPatchI,        // do not use patch info
+                    edgeZoneID[edgeI],
+                    edgeFlip[edgeI]
+                );
+            }
+            else
+            {
+                // Proper, uncoupled patch edge
+
+                findZoneFace
+                (
+                    false,      // useInternalFaces,
+                    true,       // useBoundaryFaces,
+
+                    mesh,
+                    pp,
+                    edgeI,
+
+
+                    ppFaces,    // excludeFaces,
+                    meshFaces,  // meshFaces,
+
+                    inflateFaceID[edgeI],
+                    edgePatchID[edgeI],
+                    edgeZoneID[edgeI],
+                    edgeFlip[edgeI]
+                );
             }
         }
     }
 
 
 
-    // Now hopefully every boundary edge has a side patch. Check
+    // Now hopefully every boundary edge has a edge patch. Check
     if (debug)
     {
         forAll(edgeFaces, edgeI)
         {
-            if (edgeFaces[edgeI].size() == 1 && sidePatchID[edgeI] == -1)
+            if (edgeFaces[edgeI].size() == 1 && edgePatchID[edgeI] == -1)
             {
                 const edge& e = pp.edges()[edgeI];
-                //FatalErrorIn("addPatchCellLayer::calcSidePatch(..)")
-                WarningIn("addPatchCellLayer::calcSidePatch(..)")
-                    << "Have no sidePatchID for edge " << edgeI << " points "
+                //FatalErrorIn("addPatchCellLayer::calcExtrudeInfo(..)")
+                WarningIn("addPatchCellLayer::calcExtrudeInfo(..)")
+                    << "Have no edgePatchID for edge " << edgeI << " points "
                     << pp.points()[pp.meshPoints()[e[0]]]
                     << pp.points()[pp.meshPoints()[e[1]]]
                     //<< abort(FatalError);
@@ -718,15 +866,15 @@ void Foam::addPatchCellLayer::calcSidePatch
 
 
 
-    // Now we have sidepatch see if we have patchface or edge to inflate
-    // from.
+    // Pass3: for any faces set in pass1 see if we can find a processor face
+    //        to inherit from (we only have a patch, not a patch face)
     forAll(edgeFaces, edgeI)
     {
         if
         (
             edgeFaces[edgeI].size() == 1
-         && sidePatchID[edgeI] != -1
-         && inflateFaceI[edgeI] == -1
+         && edgePatchID[edgeI] != -1
+         && inflateFaceID[edgeI] == -1
         )
         {
             // 1. Do we have a boundary face to inflate from
@@ -745,35 +893,116 @@ void Foam::addPatchCellLayer::calcSidePatch
             {
                 label faceI = meshFaces[k];
 
-                if (faceI != myFaceI)
+                if (faceI != myFaceI && !mesh.isInternalFace(faceI))
                 {
-                    if (mesh.isInternalFace(faceI))
+                    if (patches.whichPatch(faceI) == edgePatchID[edgeI])
                     {
-                        inflateEdgeI[edgeI] = meshEdgeI;
-                    }
-                    else
-                    {
-                        if (patches.whichPatch(faceI) == sidePatchID[edgeI])
-                        {
-                            setFaceProps
-                            (
-                                mesh,
-                                faceI,
+                        setFaceProps
+                        (
+                            mesh,
+                            pp,
+                            edgeI,
+                            faceI,
 
-                                sidePatchID[edgeI],
-                                sideZoneID[edgeI],
-                                sideFlip[edgeI]
-                            );
-                            inflateFaceI[edgeI] = faceI;
-                            inflateEdgeI[edgeI] = -1;
-
-                            break;
-                        }
+                            edgePatchID[edgeI],
+                            edgeZoneID[edgeI],
+                            edgeFlip[edgeI],
+                            inflateFaceID[edgeI]
+                        );
+                        break;
                     }
                 }
             }
         }
     }
+
+
+
+    // Sync all data:
+    // - edgePatchID : might be local in case of processor patch. Do not
+    //   sync for now
+    // - inflateFaceID: local. Do not sync
+    // - edgeZoneID : global. sync.
+    // - edgeFlip : global. sync.
+
+    if (Pstream::parRun())
+    {
+        const globalMeshData& gd = mesh.globalData();
+        const indirectPrimitivePatch& cpp = gd.coupledPatch();
+
+        labelList patchEdges;
+        labelList coupledEdges;
+        PackedBoolList sameEdgeOrientation;
+        PatchTools::matchEdges
+        (
+            pp,
+            cpp,
+            patchEdges,
+            coupledEdges,
+            sameEdgeOrientation
+        );
+
+        // Convert data on pp edges to data on coupled patch
+        labelList cppEdgeZoneID(cpp.nEdges(), -1);
+        boolList cppEdgeFlip(cpp.nEdges(), false);
+        forAll(coupledEdges, i)
+        {
+            label cppEdgeI = coupledEdges[i];
+            label ppEdgeI = patchEdges[i];
+
+            cppEdgeZoneID[cppEdgeI] = edgeZoneID[ppEdgeI];
+            if (sameEdgeOrientation[i])
+            {
+                cppEdgeFlip[cppEdgeI] = edgeFlip[ppEdgeI];
+            }
+            else
+            {
+                cppEdgeFlip[cppEdgeI] = !edgeFlip[ppEdgeI];
+            }
+        }
+
+        // Sync
+        const globalIndexAndTransform& git = gd.globalTransforms();
+        const mapDistribute& edgeMap = gd.globalEdgeSlavesMap();
+
+        globalMeshData::syncData
+        (
+            cppEdgeZoneID,
+            gd.globalEdgeSlaves(),
+            gd.globalEdgeTransformedSlaves(),
+            edgeMap,
+            git,
+            maxEqOp<label>(),
+            dummyTransform()
+        );
+        globalMeshData::syncData
+        (
+            cppEdgeFlip,
+            gd.globalEdgeSlaves(),
+            gd.globalEdgeTransformedSlaves(),
+            edgeMap,
+            git,
+            andEqOp<bool>(),
+            dummyTransform()
+        );
+
+        // Convert data on coupled edges to pp edges
+        forAll(coupledEdges, i)
+        {
+            label cppEdgeI = coupledEdges[i];
+            label ppEdgeI = patchEdges[i];
+
+            edgeZoneID[ppEdgeI] = cppEdgeZoneID[cppEdgeI];
+            if (sameEdgeOrientation[i])
+            {
+                edgeFlip[ppEdgeI] = cppEdgeFlip[cppEdgeI];
+            }
+            else
+            {
+                edgeFlip[ppEdgeI] = !cppEdgeFlip[cppEdgeI];
+            }
+        }
+    }
 }
 
 
@@ -783,7 +1012,12 @@ void Foam::addPatchCellLayer::setRefinement
     const labelListList& globalEdgeFaces,
     const scalarField& expansionRatio,
     const indirectPrimitivePatch& pp,
-    const labelList& sidePatchID,
+
+    const labelList& edgePatchID,
+    const labelList& edgeZoneID,
+    const boolList& edgeFlip,
+    const labelList& inflateFaceID,
+
     const labelList& exposedPatchID,
     const labelList& nFaceLayers,
     const labelList& nPointLayers,
@@ -865,7 +1099,10 @@ void Foam::addPatchCellLayer::setRefinement
                     << e.line(pp.localPoints())
                     << " which is non-manifold (has "
                     << globalEdgeFaces[edgeI].size()
-                    << " faces using it)"
+                    << " local or coupled faces using it)"
+                    << " of which "
+                    << pp.edgeFaces()[edgeI].size()
+                    << " local"
                     << abort(FatalError);
             }
         }
@@ -1741,13 +1978,25 @@ void Foam::addPatchCellLayer::setRefinement
                             ef
                         );
 
+                        // Because we walk in order of patch face and in order
+                        // of face edges so face orientation will be opposite
+                        // that of the patch edge
+                        bool zoneFlip = false;
+                        if (edgeZoneID[startEdgeI] != -1)
+                        {
+                            zoneFlip = !edgeFlip[startEdgeI];
+                        }
+
                         addSideFace
                         (
                             pp,
                             addedCells,
 
                             newFace,                // vertices of new face
-                            sidePatchID[startEdgeI],// -1 or patch for face
+                            edgePatchID[startEdgeI],// -1 or patch for face
+                            edgeZoneID[startEdgeI],
+                            zoneFlip,
+                            inflateFaceID[startEdgeI],
 
                             patchFaceI,
                             nbrFaceI,

src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -220,6 +220,9 @@ class addPatchCellLayer
 
             const face& newFace,
             const label newPatchID,
+            const label newZoneI,
+            const bool newFlip,
+            const label inflateFaceI,
 
             const label ownFaceI,
             const label nbrFaceI,
@@ -243,6 +246,39 @@ class addPatchCellLayer
             bool&
         );
 
+        //- Extract properties from mesh face in pp edge ordering
+        static void setFaceProps
+        (
+            const polyMesh& mesh,
+            const indirectPrimitivePatch& pp,
+            const label ppEdgeI,
+            const label faceI,
+
+            label& patchI,
+            label& zoneI,
+            bool& zoneFlip,
+            label& inflateFaceI
+        );
+
+        //- Find internal or boundary face to get extrude properties
+        //  from. zoneFlip consistent with ppEdge ordering
+        static void findZoneFace
+        (
+            const bool useInternalFaces,
+            const bool useBoundaryFaces,
+
+            const polyMesh& mesh,
+            const indirectPrimitivePatch& pp,
+            const label ppEdgeI,
+            const UIndirectList<label>& excludeFaces,
+            const labelList& meshFaces,
+
+            label& inflateFaceI,
+            label& patchI,
+            label& zoneI,
+            bool& zoneFlip
+        );
+
         //- Disallow default bitwise copy construct
         addPatchCellLayer(const addPatchCellLayer&);
 
@@ -303,22 +339,37 @@ public:
                 const indirectPrimitivePatch& pp
             );
 
-            //- Boundary edges get extruded into boundary faces. Determine patch
-            //  for these faces. This might be a to-be-created processor patch
-            //  (patchI >= mesh.boundaryMesh().size()) in which case the
-            //  nbrProcToPatch, patchToNbrProc give the correspondence. nPatches
-            //  is the new number of patches.
-            static void calcSidePatch
+            //- Determine extrude information per patch edge:
+            //  - zoneID, zoneFlip :
+            //      picks one of the faces that connects to
+            //      the edge. For boundary edge only looks
+            //      at boundary faces. For internal edge it looks at internal
+            //      faces only (zoneFromAnyFace = false) or at any face
+            //      (zoneFromAnyFace = true). zoneFlip is consistent with
+            //      ordering of pp edge.
+            //      Face selected gets stored in inflateFaceID
+            //  - patchID :
+            //      get patch from any boundary face connected to the
+            //      edge. The patch might be a to-be-created processor patch
+            //      (patchI >= mesh.boundaryMesh().size()) in which case the
+            //      nbrProcToPatch, patchToNbrProc give the correspondence.
+            //      nPatches is the new number of patches.
+            static void calcExtrudeInfo
             (
+                const bool zoneFromAnyFace,
+
                 const polyMesh&,
                 const globalIndex& globalFaces,
                 const labelListList& globalEdgeFaces,
                 const indirectPrimitivePatch& pp,
 
-                labelList& sidePatchID,
+                labelList& edgePatchID,     // if extruding a patch edge
                 label& nPatches,
                 Map<label>& nbrProcToPatch,
-                Map<label>& patchToNbrProc
+                Map<label>& patchToNbrProc,
+                labelList& edgeZoneID,
+                boolList& edgeFlip,
+                labelList& inflateFaceID
             );
 
             //- Play commands into polyTopoChange to create layers on top
@@ -347,7 +398,12 @@ public:
                 const labelListList& globalEdgeFaces,
                 const scalarField& expansionRatio,
                 const indirectPrimitivePatch& pp,
+
                 const labelList& sidePatchID,
+                const labelList& sideZoneID,
+                const boolList& sideFlip,
+                const labelList& inflateFaceID,
+
                 const labelList& exposedPatchID,
                 const labelList& nFaceLayers,
                 const labelList& nPointLayers,
@@ -374,7 +430,12 @@ public:
                     globalEdgeFaces,
                     scalarField(pp.nPoints(), 1.0),     // expansion ration
                     pp,
+
                     sidePatchID,
+                    labelList(pp.nEdges(), -1),         // zoneID
+                    boolList(pp.nEdges(), false),       // zoneFlip
+                    labelList(pp.nEdges(), -1),         // inflateFaceID
+
                     labelList(0),
                     labelList(pp.size(), nLayers),      // nFaceLayers
                     labelList(pp.nPoints(), nLayers),   // nPointLayers

src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -795,7 +795,7 @@ Foam::label Foam::hexRef8::findLevel
 
 
 // Gets cell level such that the face has four points <= level.
-Foam::label Foam::hexRef8::getAnchorLevel(const label faceI) const
+Foam::label Foam::hexRef8::faceLevel(const label faceI) const
 {
     const face& f = mesh_.faces()[faceI];
 
@@ -3519,7 +3519,7 @@ Foam::labelListList Foam::hexRef8::setRefinement
 
     for (label faceI = 0; faceI < mesh_.nFaces(); faceI++)
     {
-        faceAnchorLevel[faceI] = getAnchorLevel(faceI);
+        faceAnchorLevel[faceI] = faceLevel(faceI);
     }
 
     // -1  : no need to split face

src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -411,7 +411,7 @@ public:
         // Refinement
 
             //- Gets level such that the face has four points <= level.
-            label getAnchorLevel(const label faceI) const;
+            label faceLevel(const label faceI) const;
 
             //- Given valid mesh and current cell level and proposed
             //  cells to refine calculate any clashes (due to 2:1) and return

src/mesh/autoMesh/Make/files

@@ -2,7 +2,6 @@ autoHexMesh             = autoHexMesh
 autoHexMeshDriver       = $(autoHexMesh)/autoHexMeshDriver
 
 $(autoHexMeshDriver)/autoLayerDriver.C
-/* $(autoHexMeshDriver)/autoLayerDriverShrink.C */
 $(autoHexMeshDriver)/autoSnapDriver.C
 $(autoHexMeshDriver)/autoSnapDriverFeature.C
 $(autoHexMeshDriver)/autoRefineDriver.C
@@ -10,7 +9,6 @@ $(autoHexMeshDriver)/autoRefineDriver.C
 $(autoHexMeshDriver)/layerParameters/layerParameters.C
 $(autoHexMeshDriver)/refinementParameters/refinementParameters.C
 $(autoHexMeshDriver)/snapParameters/snapParameters.C
-$(autoHexMeshDriver)/pointData/pointData.C
 
 $(autoHexMesh)/meshRefinement/meshRefinementBaffles.C
 $(autoHexMesh)/meshRefinement/meshRefinement.C
@@ -30,7 +28,10 @@ meshMover = $(autoHexMesh)/externalDisplacementMeshMover
 $(meshMover)/displacementMeshMoverMotionSolver.C
 $(meshMover)/externalDisplacementMeshMover.C
 $(meshMover)/medialAxisMeshMover.C
+$(meshMover)/displacementMotionSolverMeshMover.C
+/* $(meshMover)/pointSmoothingMeshMover.C */
 $(meshMover)/zeroFixedValue/zeroFixedValuePointPatchFields.C
+$(meshMover)/fieldSmoother/fieldSmoother.C
 
 LIB = $(FOAM_LIBBIN)/libautoMesh
 

src/mesh/autoMesh/Make/options

@@ -18,4 +18,5 @@ LIB_LIBS = \
     -ledgeMesh \
     -lsurfMesh \
     -ltriSurface \
+    -lfvMotionSolvers \
     -ldistributed

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -47,7 +47,6 @@ class removePoints;
 class pointSet;
 class motionSmoother;
 class addPatchCellLayer;
-class pointData;
 class faceSet;
 class layerParameters;
 
@@ -108,7 +107,6 @@ private:
         const labelList globalToSlavePatch_;
 
 
-
     // Private Member Functions
 
         // Layers
@@ -231,14 +229,17 @@ private:
                     List<extrudeMode>& extrudeStatus
                 ) const;
 
-                //- See what patches boundaryedges should be extruded into
+                //- See what zones and patches edges should be extruded into
                 void determineSidePatches
                 (
                     const globalIndex& globalFaces,
                     const labelListList& edgeGlobalFaces,
                     const indirectPrimitivePatch& pp,
 
-                    labelList& sidePatchID
+                    labelList& edgePatchID,
+                    labelList& edgeZoneID,
+                    boolList& edgeFlip,
+                    labelList& inflateFaceID
                 );
 
                 //- Calculate pointwise wanted and minimum thickness.
@@ -370,6 +371,14 @@ private:
                     const List<extrudeMode>& extrudeStatus
                 );
 
+                //- After adding to mesh get the new baffles
+                static List<labelPair> getBafflesOnAddedMesh
+                (
+                    const polyMesh& mesh,
+                    const labelList& newToOldFaces,
+                    const List<labelPair>& baffles
+                );
+
                 //- Collect layer faces and layer cells into bools
                 //  for ease of handling
                 static void getLayerCellsFaces
@@ -393,6 +402,15 @@ private:
                 ) const;
 
                 //- Write cellSet,faceSet for layers
+                bool writeLayerSets
+                (
+                    const fvMesh& mesh,
+                    const labelList& cellNLayers,
+                    const scalarField& faceRealThickness
+                ) const;
+
+                //- Write volFields,cellSet,faceSet for layers depending
+                //  on write level
                 bool writeLayerData
                 (
                     const fvMesh& mesh,
@@ -464,13 +482,6 @@ private:
                     pointVectorField& normals
                 ) const;
 
-                bool isMaxEdge
-                (
-                    const List<pointData>&,
-                    const label edgeI,
-                    const scalar minCos
-                ) const;
-
                 //- Stop layer growth where mesh wraps around edge with a
                 //  large feature angle
                 void handleFeatureAngleLayerTerminations
@@ -595,6 +606,7 @@ public:
                 const dictionary& shrinkDict,
                 const dictionary& motionDict,
                 const layerParameters& layerParams,
+                const bool mergePatchFaces,         // merging patch faces
                 const bool preBalance,              // balance before adding?
                 decompositionMethod& decomposer,
                 fvMeshDistribute& distributor

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,6 +37,8 @@ License
 #include "unitConversion.H"
 #include "snapParameters.H"
 #include "localPointRegion.H"
+#include "IOmanip.H"
+#include "labelVector.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -94,7 +96,7 @@ Foam::label Foam::autoRefineDriver::featureEdgeRefine
             (
                 meshRefiner_.refineCandidates
                 (
-                    refineParams.keepPoints(),
+                    refineParams.locationsInMesh(),
                     refineParams.curvature(),
                     refineParams.planarAngle(),
 
@@ -207,7 +209,7 @@ Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
         (
             meshRefiner_.refineCandidates
             (
-                refineParams.keepPoints(),
+                refineParams.locationsInMesh(),
                 refineParams.curvature(),
                 refineParams.planarAngle(),
 
@@ -341,7 +343,7 @@ Foam::label Foam::autoRefineDriver::gapOnlyRefine
         (
             meshRefiner_.refineCandidates
             (
-                refineParams.keepPoints(),
+                refineParams.locationsInMesh(),
                 refineParams.curvature(),
                 refineParams.planarAngle(),
 
@@ -669,6 +671,348 @@ Foam::label Foam::autoRefineDriver::danglingCellRefine
 }
 
 
+// Detect cells with opposing intersected faces of differing refinement
+// level and refine them.
+Foam::label Foam::autoRefineDriver::refinementInterfaceRefine
+(
+    const refinementParameters& refineParams,
+    const label maxIter
+)
+{
+    const fvMesh& mesh = meshRefiner_.mesh();
+
+    label iter = 0;
+
+    if (refineParams.interfaceRefine())
+    {
+        for (;iter < maxIter; iter++)
+        {
+            Info<< nl
+                << "Refinement transition refinement iteration " << iter << nl
+                << "--------------------------------------------" << nl
+                << endl;
+
+            const labelList& surfaceIndex = meshRefiner_.surfaceIndex();
+            const hexRef8& cutter = meshRefiner_.meshCutter();
+            const vectorField& fA = mesh.faceAreas();
+            const labelList& faceOwner = mesh.faceOwner();
+
+
+            // Determine cells to refine
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+            const cellList& cells = mesh.cells();
+
+            labelList candidateCells;
+            {
+                // Pass1: pick up cells with differing face level
+
+                cellSet transitionCells
+                (
+                    mesh,
+                    "transitionCells",
+                    cells.size()/100
+                );
+
+                forAll(cells, cellI)
+                {
+                    const cell& cFaces = cells[cellI];
+                    label cLevel = cutter.cellLevel()[cellI];
+
+                    forAll(cFaces, cFaceI)
+                    {
+                        label faceI = cFaces[cFaceI];
+
+                        if (surfaceIndex[faceI] != -1)
+                        {
+                            label fLevel = cutter.faceLevel(faceI);
+                            if (fLevel != cLevel)
+                            {
+                                transitionCells.insert(cellI);
+                            }
+                        }
+                    }
+                }
+
+
+                cellSet candidateCellSet
+                (
+                    mesh,
+                    "candidateCells",
+                    cells.size()/1000
+                );
+
+                // Pass2: check for oppositeness
+
+                //forAllConstIter(cellSet, transitionCells, iter)
+                //{
+                //    label cellI = iter.key();
+                //    const cell& cFaces = cells[cellI];
+                //    const point& cc = cellCentres[cellI];
+                //    const scalar rCVol = pow(cellVolumes[cellI], -5.0/3.0);
+                //
+                //    // Determine principal axes of cell
+                //    symmTensor R(symmTensor::zero);
+                //
+                //    forAll(cFaces, i)
+                //    {
+                //        label faceI = cFaces[i];
+                //
+                //        const point& fc = faceCentres[faceI];
+                //
+                //        // Calculate face-pyramid volume
+                //        scalar pyrVol = 1.0/3.0 * fA[faceI] & (fc-cc);
+                //
+                //        if (faceOwner[faceI] != cellI)
+                //        {
+                //            pyrVol = -pyrVol;
+                //        }
+                //
+                //        // Calculate face-pyramid centre
+                //        vector pc = (3.0/4.0)*fc + (1.0/4.0)*cc;
+                //
+                //        R += pyrVol*sqr(pc-cc)*rCVol;
+                //    }
+                //
+                //    //- MEJ: Problem: truncation errors cause complex evs
+                //    vector lambdas(eigenValues(R));
+                //    const tensor axes(eigenVectors(R, lambdas));
+                //
+                //
+                //    // Check if this cell has
+                //    // - opposing sides intersected
+                //    // - which are of different refinement level
+                //    // - plus the inbetween face
+                //
+                //    labelVector plusFaceLevel(labelVector(-1, -1, -1));
+                //    labelVector minFaceLevel(labelVector(-1, -1, -1));
+                //
+                //    forAll(cFaces, cFaceI)
+                //    {
+                //        label faceI = cFaces[cFaceI];
+                //
+                //        if (surfaceIndex[faceI] != -1)
+                //        {
+                //            label fLevel = cutter.faceLevel(faceI);
+                //
+                //            // Get outwards pointing normal
+                //            vector n = fA[faceI]/mag(fA[faceI]);
+                //            if (faceOwner[faceI] != cellI)
+                //            {
+                //                n = -n;
+                //            }
+                //
+                //            // What is major direction and sign
+                //            direction cmpt = vector::X;
+                //            scalar maxComp = (n&axes.x());
+                //
+                //            scalar yComp = (n&axes.y());
+                //            scalar zComp = (n&axes.z());
+                //
+                //            if (mag(yComp) > mag(maxComp))
+                //            {
+                //                maxComp = yComp;
+                //                cmpt = vector::Y;
+                //            }
+                //
+                //            if (mag(zComp) > mag(maxComp))
+                //            {
+                //                maxComp = zComp;
+                //                cmpt = vector::Z;
+                //            }
+                //
+                //            if (maxComp > 0)
+                //            {
+                //                plusFaceLevel[cmpt] = max
+                //                (
+                //                    plusFaceLevel[cmpt],
+                //                    fLevel
+                //                );
+                //            }
+                //            else
+                //            {
+                //                minFaceLevel[cmpt] = max
+                //                (
+                //                    minFaceLevel[cmpt],
+                //                    fLevel
+                //                );
+                //            }
+                //        }
+                //    }
+                //
+                //    // Check if we picked up any opposite differing level
+                //    for (direction dir = 0; dir < vector::nComponents; dir++)
+                //    {
+                //        if
+                //        (
+                //            plusFaceLevel[dir] != -1
+                //         && minFaceLevel[dir] != -1
+                //         && plusFaceLevel[dir] != minFaceLevel[dir]
+                //        )
+                //        {
+                //            candidateCellSet.insert(cellI);
+                //        }
+                //    }
+                //}
+
+                const scalar oppositeCos = Foam::cos(Foam::degToRad(135));
+
+                forAllConstIter(cellSet, transitionCells, iter)
+                {
+                    label cellI = iter.key();
+                    const cell& cFaces = cells[cellI];
+                    label cLevel = cutter.cellLevel()[cellI];
+
+                    // Detect opposite intersection
+                    bool foundOpposite = false;
+
+                    forAll(cFaces, cFaceI)
+                    {
+                        label faceI = cFaces[cFaceI];
+
+                        if
+                        (
+                            surfaceIndex[faceI] != -1
+                         && cutter.faceLevel(faceI) > cLevel
+                        )
+                        {
+                            // Get outwards pointing normal
+                            vector n = fA[faceI]/mag(fA[faceI]);
+                            if (faceOwner[faceI] != cellI)
+                            {
+                                n = -n;
+                            }
+
+                            // Check for any opposite intersection
+                            forAll(cFaces, cFaceI2)
+                            {
+                                label face2I = cFaces[cFaceI2];
+
+                                if
+                                (
+                                    face2I != faceI
+                                 && surfaceIndex[face2I] != -1
+                                )
+                                {
+                                    // Get outwards pointing normal
+                                    vector n2 = fA[face2I]/mag(fA[face2I]);
+                                    if (faceOwner[face2I] != cellI)
+                                    {
+                                        n2 = -n2;
+                                    }
+
+
+                                    if ((n&n2) < oppositeCos)
+                                    {
+                                        foundOpposite = true;
+                                        break;
+                                    }
+                                }
+                            }
+
+                            if (foundOpposite)
+                            {
+                                break;
+                            }
+                        }
+                    }
+
+
+                    if (foundOpposite)
+                    {
+                        candidateCellSet.insert(cellI);
+                    }
+                }
+
+                if (debug&meshRefinement::MESH)
+                {
+                    Pout<< "Dumping " << candidateCellSet.size()
+                        << " cells to cellSet candidateCellSet." << endl;
+                    candidateCellSet.instance() = meshRefiner_.timeName();
+                    candidateCellSet.write();
+                }
+                candidateCells = candidateCellSet.toc();
+            }
+
+
+
+            labelList cellsToRefine
+            (
+                meshRefiner_.meshCutter().consistentRefinement
+                (
+                    candidateCells,
+                    true
+                )
+            );
+            Info<< "Determined cells to refine in = "
+                << mesh.time().cpuTimeIncrement() << " s" << endl;
+
+
+            label nCellsToRefine = cellsToRefine.size();
+            reduce(nCellsToRefine, sumOp<label>());
+
+            Info<< "Selected for refinement : " << nCellsToRefine
+                << " cells (out of " << mesh.globalData().nTotalCells()
+                << ')' << endl;
+
+            // Stop when no cells to refine. After a few iterations check if too
+            // few cells
+            if
+            (
+                nCellsToRefine == 0
+             || (
+                    iter >= 1
+                 && nCellsToRefine <= refineParams.minRefineCells()
+                )
+            )
+            {
+                Info<< "Stopping refining since too few cells selected."
+                    << nl << endl;
+                break;
+            }
+
+
+            if (debug)
+            {
+                const_cast<Time&>(mesh.time())++;
+            }
+
+
+            if
+            (
+                returnReduce
+                (
+                    (mesh.nCells() >= refineParams.maxLocalCells()),
+                    orOp<bool>()
+                )
+            )
+            {
+                meshRefiner_.balanceAndRefine
+                (
+                    "interface cell refinement iteration " + name(iter),
+                    decomposer_,
+                    distributor_,
+                    cellsToRefine,
+                    refineParams.maxLoadUnbalance()
+                );
+            }
+            else
+            {
+                meshRefiner_.refineAndBalance
+                (
+                    "interface cell refinement iteration " + name(iter),
+                    decomposer_,
+                    distributor_,
+                    cellsToRefine,
+                    refineParams.maxLoadUnbalance()
+                );
+            }
+        }
+    }
+    return iter;
+}
+
+
 void Foam::autoRefineDriver::removeInsideCells
 (
     const refinementParameters& refineParams,
@@ -692,13 +1036,15 @@ void Foam::autoRefineDriver::removeInsideCells
         nBufferLayers,                  // nBufferLayers
         globalToMasterPatch_,
         globalToSlavePatch_,
-        refineParams.keepPoints()[0]
+        refineParams.locationsInMesh(),
+        refineParams.zonesInMesh(),
+        refineParams.locationsOutsideMesh()
     );
 
     if (debug&meshRefinement::MESH)
     {
         Pout<< "Writing subsetted mesh to time "
-            << meshRefiner_.timeName() << '.' << endl;
+            << meshRefiner_.timeName() << endl;
         meshRefiner_.write
         (
             meshRefinement::debugType(debug),
@@ -753,7 +1099,7 @@ Foam::label Foam::autoRefineDriver::shellRefine
         (
             meshRefiner_.refineCandidates
             (
-                refineParams.keepPoints(),
+                refineParams.locationsInMesh(),
                 refineParams.curvature(),
                 refineParams.planarAngle(),
 
@@ -913,22 +1259,40 @@ void Foam::autoRefineDriver::baffleAndSplitMesh
         false,                          // perpendicular edge connected cells
         scalarField(0),                 // per region perpendicular angle
 
-        // Free standing baffles
-        !handleSnapProblems,            // merge free standing baffles?
-        refineParams.planarAngle(),
-
         motionDict,
         const_cast<Time&>(mesh.time()),
         globalToMasterPatch_,
         globalToSlavePatch_,
-        refineParams.keepPoints()[0]
+        refineParams.locationsInMesh(),
+        refineParams.zonesInMesh(),
+        refineParams.locationsOutsideMesh()
     );
+
+
+    if (!handleSnapProblems) // merge free standing baffles?
+    {
+        meshRefiner_.mergeFreeStandingBaffles
+        (
+            snapParams,
+            refineParams.useTopologicalSnapDetection(),
+            false,                  // perpendicular edge connected cells
+            scalarField(0),         // per region perpendicular angle
+            refineParams.planarAngle(),
+            motionDict,
+            const_cast<Time&>(mesh.time()),
+            globalToMasterPatch_,
+            globalToSlavePatch_,
+            refineParams.locationsInMesh(),
+            refineParams.locationsOutsideMesh()
+        );
+    }
 }
 
 
 void Foam::autoRefineDriver::zonify
 (
-    const refinementParameters& refineParams
+    const refinementParameters& refineParams,
+    wordPairHashTable& zonesToFaceZone
 )
 {
     // Mesh is at its finest. Do zoning
@@ -940,7 +1304,11 @@ void Foam::autoRefineDriver::zonify
     const labelList namedSurfaces =
         surfaceZonesInfo::getNamedSurfaces(meshRefiner_.surfaces().surfZones());
 
-    if (namedSurfaces.size())
+    if
+    (
+        namedSurfaces.size()
+     || refineParams.zonesInMesh().size()
+    )
     {
         Info<< nl
             << "Introducing zones for interfaces" << nl
@@ -956,14 +1324,16 @@ void Foam::autoRefineDriver::zonify
 
         meshRefiner_.zonify
         (
-            refineParams.keepPoints()[0],
-            refineParams.allowFreeStandingZoneFaces()
+            refineParams.allowFreeStandingZoneFaces(),
+            refineParams.locationsInMesh(),
+            refineParams.zonesInMesh(),
+            zonesToFaceZone
         );
 
         if (debug&meshRefinement::MESH)
         {
             Pout<< "Writing zoned mesh to time "
-                << meshRefiner_.timeName() << '.' << endl;
+                << meshRefiner_.timeName() << endl;
             meshRefiner_.write
             (
                 meshRefinement::debugType(debug),
@@ -1015,15 +1385,29 @@ void Foam::autoRefineDriver::splitAndMergeBaffles
         handleSnapProblems,                 // remove perp edge connected cells
         perpAngle,                          // perp angle
 
-        // Free standing baffles
-        true,                               // merge free standing baffles?
-        refineParams.planarAngle(),         // planar angle
-
         motionDict,
         const_cast<Time&>(mesh.time()),
         globalToMasterPatch_,
         globalToSlavePatch_,
-        refineParams.keepPoints()[0]
+        refineParams.locationsInMesh(),
+        refineParams.zonesInMesh(),
+        refineParams.locationsOutsideMesh()
+    );
+
+    // Merge free-standing baffles always
+    meshRefiner_.mergeFreeStandingBaffles
+    (
+        snapParams,
+        refineParams.useTopologicalSnapDetection(),
+        handleSnapProblems,
+        perpAngle,
+        refineParams.planarAngle(),
+        motionDict,
+        const_cast<Time&>(mesh.time()),
+        globalToMasterPatch_,
+        globalToSlavePatch_,
+        refineParams.locationsInMesh(),
+        refineParams.locationsOutsideMesh()
     );
 
     if (debug)
@@ -1048,7 +1432,7 @@ void Foam::autoRefineDriver::splitAndMergeBaffles
         // Actually merge baffles. Note: not exactly parallellized. Should
         // convert baffle faces into processor faces if they resulted
         // from them.
-        meshRefiner_.mergeBaffles(couples);
+        meshRefiner_.mergeBaffles(couples, Map<label>(0));
 
         if (debug)
         {
@@ -1061,7 +1445,8 @@ void Foam::autoRefineDriver::splitAndMergeBaffles
         (
             globalToMasterPatch_,
             globalToSlavePatch_,
-            refineParams.keepPoints()[0]
+            refineParams.locationsInMesh(),
+            refineParams.locationsOutsideMesh()
         );
 
         if (debug)
@@ -1077,7 +1462,7 @@ void Foam::autoRefineDriver::splitAndMergeBaffles
     if (debug&meshRefinement::MESH)
     {
         Pout<< "Writing handleProblemCells mesh to time "
-            << meshRefiner_.timeName() << '.' << endl;
+            << meshRefiner_.timeName() << endl;
         meshRefiner_.write
         (
             meshRefinement::debugType(debug),
@@ -1092,8 +1477,83 @@ void Foam::autoRefineDriver::splitAndMergeBaffles
 }
 
 
+void Foam::autoRefineDriver::addFaceZones
+(
+    meshRefinement& meshRefiner,
+    const refinementParameters& refineParams,
+    const HashTable<Pair<word> >& faceZoneToPatches
+)
+{
+    if (faceZoneToPatches.size())
+    {
+        Info<< nl
+            << "Adding patches for face zones" << nl
+            << "-----------------------------" << nl
+            << endl;
+
+        Info<< setf(ios_base::left)
+            << setw(6) << "Patch"
+            << setw(20) << "Type"
+            << setw(30) << "Name"
+            << setw(30) << "FaceZone"
+            << setw(10) << "FaceType"
+            << nl
+            << setw(6) << "-----"
+            << setw(20) << "----"
+            << setw(30) << "----"
+            << setw(30) << "--------"
+            << setw(10) << "--------"
+            << endl;
+
+        const polyMesh& mesh = meshRefiner.mesh();
+
+        // Add patches for added inter-region faceZones
+        forAllConstIter(HashTable<Pair<word> >, faceZoneToPatches, iter)
+        {
+            const word& fzName = iter.key();
+            const Pair<word>& patchNames = iter();
+
+            // Get any user-defined faceZone data
+            surfaceZonesInfo::faceZoneType fzType;
+            dictionary patchInfo = refineParams.getZoneInfo(fzName, fzType);
+
+            const word& masterName = fzName;
+            //const word slaveName = fzName + "_slave";
+            //const word slaveName = czNames.second()+"_to_"+czNames.first();
+            const word& slaveName = patchNames.second();
+
+            label mpI = meshRefiner.addMeshedPatch(masterName, patchInfo);
+
+            Info<< setf(ios_base::left)
+                << setw(6) << mpI
+                << setw(20) << mesh.boundaryMesh()[mpI].type()
+                << setw(30) << masterName
+                << setw(30) << fzName
+                << setw(10) << surfaceZonesInfo::faceZoneTypeNames[fzType]
+                << nl;
+
+
+            label slI = meshRefiner.addMeshedPatch(slaveName, patchInfo);
+
+            Info<< setf(ios_base::left)
+                << setw(6) << slI
+                << setw(20) << mesh.boundaryMesh()[slI].type()
+                << setw(30) << slaveName
+                << setw(30) << fzName
+                << setw(10) << surfaceZonesInfo::faceZoneTypeNames[fzType]
+                << nl;
+
+            meshRefiner.addFaceZone(fzName, masterName, slaveName, fzType);
+        }
+
+        Info<< endl;
+    }
+}
+
+
 void Foam::autoRefineDriver::mergePatchFaces
 (
+    const bool geometricMerge,
     const refinementParameters& refineParams,
     const dictionary& motionDict
 )
@@ -1105,14 +1565,28 @@ void Foam::autoRefineDriver::mergePatchFaces
 
     const fvMesh& mesh = meshRefiner_.mesh();
 
-    meshRefiner_.mergePatchFacesUndo
-    (
-        Foam::cos(degToRad(45.0)),
-        Foam::cos(degToRad(45.0)),
-        meshRefiner_.meshedPatches(),
-        motionDict,
-        labelList(mesh.nFaces(), -1)
-    );
+    if (geometricMerge)
+    {
+        meshRefiner_.mergePatchFacesUndo
+        (
+            Foam::cos(degToRad(45.0)),
+            Foam::cos(degToRad(45.0)),
+            meshRefiner_.meshedPatches(),
+            motionDict,
+            labelList(mesh.nFaces(), -1)
+        );
+    }
+    else
+    {
+        // Still merge refined boundary faces if all four are on same patch
+        meshRefiner_.mergePatchFaces
+        (
+            Foam::cos(degToRad(45.0)),
+            Foam::cos(degToRad(45.0)),
+            4,          // only merge faces split into 4
+            meshRefiner_.meshedPatches()
+        );
+    }
 
     if (debug)
     {
@@ -1134,6 +1608,7 @@ void Foam::autoRefineDriver::doRefine
     const refinementParameters& refineParams,
     const snapParameters& snapParams,
     const bool prepareForSnapping,
+    const bool doMergePatchFaces,
     const dictionary& motionDict
 )
 {
@@ -1145,7 +1620,7 @@ void Foam::autoRefineDriver::doRefine
     const fvMesh& mesh = meshRefiner_.mesh();
 
     // Check that all the keep points are inside the mesh.
-    refineParams.findCells(mesh);
+    refineParams.findCells(true, mesh, refineParams.locationsInMesh());
 
     // Refine around feature edges
     featureEdgeRefine
@@ -1196,6 +1671,13 @@ void Foam::autoRefineDriver::doRefine
         100     // maxIter
     );
 
+    // Refine any cells with differing refinement level on either side
+    refinementInterfaceRefine
+    (
+        refineParams,
+        10      // maxIter
+    );
+
     // Introduce baffles at surface intersections. Remove sections unreachable
     // from keepPoint.
     baffleAndSplitMesh
@@ -1206,8 +1688,31 @@ void Foam::autoRefineDriver::doRefine
         motionDict
     );
 
-    // Mesh is at its finest. Do optional zoning.
-    zonify(refineParams);
+    // Mesh is at its finest. Do optional zoning (cellZones and faceZones)
+    wordPairHashTable zonesToFaceZone;
+    zonify(refineParams, zonesToFaceZone);
+
+    // Create pairs of patches for faceZones
+    {
+        HashTable<Pair<word> > faceZoneToPatches(zonesToFaceZone.size());
+
+        //    Note: zonesToFaceZone contains the same data on different
+        //          processors but in different order. We could sort the
+        //          contents but instead just loop in sortedToc order.
+        List<Pair<word> > czs(zonesToFaceZone.sortedToc());
+
+        forAll(czs, i)
+        {
+            const Pair<word>& czNames = czs[i];
+            const word& fzName = zonesToFaceZone[czNames];
+
+            const word& masterName = fzName;
+            const word slaveName = czNames.second() + "_to_" + czNames.first();
+            Pair<word> patches(masterName, slaveName);
+            faceZoneToPatches.insert(fzName, patches);
+        }
+        addFaceZones(meshRefiner_, refineParams, faceZoneToPatches);
+    }
 
     // Pull baffles apart
     splitAndMergeBaffles
@@ -1221,7 +1726,7 @@ void Foam::autoRefineDriver::doRefine
     // Do something about cells with refined faces on the boundary
     if (prepareForSnapping)
     {
-        mergePatchFaces(refineParams, motionDict);
+        mergePatchFaces(doMergePatchFaces, refineParams, motionDict);
     }
 
 
@@ -1232,28 +1737,17 @@ void Foam::autoRefineDriver::doRefine
             << "---------------------" << nl
             << endl;
 
-        //if (debug)
-        //{
-        //    const_cast<Time&>(mesh.time())++;
-        //}
-
         // Do final balancing. Keep zoned faces on one processor since the
         // snap phase will convert them to baffles and this only works for
         // internal faces.
         meshRefiner_.balance
         (
-            true,
-            false,
-            scalarField(mesh.nCells(), 1), // dummy weights
+            true,                           // keepZoneFaces
+            false,                          // keepBaffles
+            scalarField(mesh.nCells(), 1),  // cellWeights
             decomposer_,
             distributor_
         );
-
-
-        if (debug)
-        {
-            meshRefiner_.checkZoneFaces();
-        }
     }
 }
 

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -34,7 +34,8 @@ SourceFiles
 #ifndef autoRefineDriver_H
 #define autoRefineDriver_H
 
-#include "treeBoundBox.H"
+#include "wordPairHashTable.H"
+#include "labelList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -105,6 +106,13 @@ class autoRefineDriver
             const label maxIter
         );
 
+        //- Refine cells with opposite faces with differing refinement level
+        label refinementInterfaceRefine
+        (
+            const refinementParameters& refineParams,
+            const label maxIter
+        );
+
         //- Remove all cells within intersected region
         void removeInsideCells
         (
@@ -129,7 +137,11 @@ class autoRefineDriver
         );
 
         //- Add zones
-        void zonify(const refinementParameters& refineParams);
+        void zonify
+        (
+            const refinementParameters& refineParams,
+            wordPairHashTable& zonesToFaceZone
+        );
 
         void splitAndMergeBaffles
         (
@@ -142,11 +154,11 @@ class autoRefineDriver
         //- Merge refined boundary faces (from exposing coarser cell)
         void mergePatchFaces
         (
+            const bool geometricMerge,
             const refinementParameters& refineParams,
             const dictionary& motionDict
         );
 
-
         //- Disallow default bitwise copy construct
         autoRefineDriver(const autoRefineDriver&);
 
@@ -182,8 +194,18 @@ public:
             const refinementParameters& refineParams,
             const snapParameters& snapParams,
             const bool prepareForSnapping,
+            const bool mergePatchFaces,
             const dictionary& motionDict
         );
+
+        //- Helper: add faceZones and patches
+        static void addFaceZones
+        (
+            meshRefinement& meshRefiner,
+            const refinementParameters& refineParams,
+            const HashTable<Pair<word> >& faceZoneToPatches
+        );
+
 };
 
 

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,6 +37,7 @@ SourceFiles
 #define autoSnapDriver_H
 
 #include "meshRefinement.H"
+#include "DynamicField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -45,6 +46,7 @@ namespace Foam
 
 // Forward declaration of classes
 class motionSmoother;
+class refinementParameters;
 class snapParameters;
 class pointConstraint;
 
@@ -79,26 +81,34 @@ class autoSnapDriver
                 PackedBoolList&
             );
 
+            //- Calculate displacement (over all mesh points) to move points
+            //  to average of connected cell centres
+            static tmp<pointField> smoothInternalDisplacement
+            (
+                const meshRefinement& meshRefiner,
+                const motionSmoother&
+            );
+
             //- Calculate displacement per patch point to smooth out patch.
             //  Quite complicated in determining which points to move where.
-            static pointField smoothPatchDisplacement
+            static tmp<pointField> smoothPatchDisplacement
             (
                 const motionSmoother&,
                 const List<labelPair>&
             );
 
-            //static tmp<pointField> avg
-            //(
-            //    const indirectPrimitivePatch&,
-            //    const pointField&
-            //);
+            static tmp<pointField> avg
+            (
+                const indirectPrimitivePatch&,
+                const pointField&
+            );
 
             //- Calculate displacement per patch point. Wip.
-            //static tmp<pointField> smoothLambdaMuPatchDisplacement
-            //(
-            //    const motionSmoother& meshMover,
-            //    const List<labelPair>& baffles
-            //);
+            static pointField smoothLambdaMuPatchDisplacement
+            (
+                const motionSmoother& meshMover,
+                const List<labelPair>& baffles
+            );
 
 
             //- Check that face zones are synced
@@ -136,6 +146,48 @@ class autoSnapDriver
                 DynamicList<labelPair>& splits
             ) const;
 
+            //- Get per face -1 or label of opposite face if on internal/baffle
+            //  faceZone
+            labelList getInternalOrBaffleDuplicateFace() const;
+
+            //- Get points both on patch and facezone.
+            static labelList getZoneSurfacePoints
+            (
+                const fvMesh& mesh,
+                const indirectPrimitivePatch&,
+                const word& zoneName
+            );
+
+            //- Get points both on patch and facezone.
+            template<class FaceList>
+            static labelList getFacePoints
+            (
+                const indirectPrimitivePatch& pp,
+                const FaceList& faces
+            );
+
+            //- Per patch point calculate point on nearest surface.
+            //  Return displacement of patch points.
+            static void calcNearestSurface
+            (
+                const refinementSurfaces& surfaces,
+
+                const labelList& surfacesToTest,
+                const labelListList& regionsToTest,
+
+                const pointField& localPoints,
+                const labelList& zonePointIndices,
+
+                scalarField& minSnapDist,
+                labelList& snapSurf,
+                vectorField& patchDisp,
+
+                // Optional: nearest point, normal
+                pointField& nearestPoint,
+                vectorField& nearestNormal
+            );
+
+
             // Feature line snapping
 
                 //- Is point on two feature edges that make a largish angle?
@@ -177,8 +229,8 @@ class autoSnapDriver
                     const scalarField& faceSnapDist,
                     vectorField& faceDisp,
                     vectorField& faceSurfaceNormal,
-                    labelList& faceSurfaceRegion,
-                    vectorField& faceRotation
+                    labelList& faceSurfaceRegion
+                    //vectorField& faceRotation
                 ) const;
                 void calcNearestFacePointProperties
                 (
@@ -253,6 +305,86 @@ class autoSnapDriver
                     const label faceI
                 ) const;
 
+                scalar pyrVol
+                (
+                    const indirectPrimitivePatch& pp,
+                    const vectorField& featureAttraction,
+                    const face& localF,
+                    const point& cc
+                ) const;
+                void facePoints
+                (
+                    const indirectPrimitivePatch& pp,
+                    const vectorField& featureAttraction,
+                    const vectorField& nearestAttraction,
+                    const face& f,
+                    DynamicField<point>& points
+                ) const;
+                scalar pyrVol
+                (
+                    const indirectPrimitivePatch& pp,
+                    const vectorField& featureAttraction,
+                    const vectorField& nearestAttraction,
+                    const face& localF,
+                    const point& cc
+                ) const;
+                Tuple2<point, vector> centreAndNormal
+                (
+                    const indirectPrimitivePatch& pp,
+                    const vectorField& featureAttraction,
+                    const vectorField& nearestAttraction,
+                    const face& localF
+                ) const;
+                bool isSplitAlignedWithFeature
+                (
+                    const scalar featureCos,
+                    const point& newPt0,
+                    const pointConstraint& pc0,
+                    const point& newPt1,
+                    const pointConstraint& pc1
+                ) const;
+                bool isConcave
+                (
+                    const point& c0,
+                    const vector& area0,
+                    const point& c1,
+                    const vector& area1,
+                    const scalar concaveCos
+                ) const;
+                labelPair findDiagonalAttraction
+                (
+                    const scalar featureCos,
+                    const scalar concaveCos,
+                    const scalar minAreaFraction,
+                    const indirectPrimitivePatch& pp,
+                    const vectorField& patchAttraction,
+                    const List<pointConstraint>& patchConstraints,
+                    const vectorField& nearestAttraction,
+                    const vectorField& nearestNormal,
+                    const label faceI,
+
+                    DynamicField<point>& points0,
+                    DynamicField<point>& points1
+                ) const;
+
+                //- Do all logic on whether to add face cut to diagonal
+                //  attraction
+                void splitDiagonals
+                (
+                    const scalar featureCos,
+                    const scalar concaveCos,
+                    const scalar minAreaFraction,
+
+                    const indirectPrimitivePatch& pp,
+                    const vectorField& nearestAttraction,
+                    const vectorField& nearestNormal,
+
+                    vectorField& patchAttraction,
+                    List<pointConstraint>& patchConstraints,
+                    DynamicList<label>& splitFaces,
+                    DynamicList<labelPair>& splits
+                ) const;
+
                 //- Avoid attraction across face diagonal since would
                 //  cause face squeeze
                 void avoidDiagonalAttraction
@@ -264,6 +396,14 @@ class autoSnapDriver
                     List<pointConstraint>& patchConstraints
                 ) const;
 
+                //- Write some stats about constraints
+                void writeStats
+                (
+                    const indirectPrimitivePatch& pp,
+                    const PackedBoolList& isMasterPoint,
+                    const List<pointConstraint>& patchConstraints
+                ) const;
+
                 //- Return hit if on multiple points
                 pointIndexHit findMultiPatchPoint
                 (
@@ -320,7 +460,6 @@ class autoSnapDriver
                 void featureAttractionUsingReconstruction
                 (
                     const label iter,
-                    const bool avoidSnapProblems,
                     const scalar featureCos,
                     const indirectPrimitivePatch& pp,
                     const scalarField& snapDist,
@@ -366,6 +505,7 @@ class autoSnapDriver
                 void determineBaffleFeatures
                 (
                     const label iter,
+                    const bool baffleFeaturePoints,
                     const scalar featureCos,
 
                     const indirectPrimitivePatch& pp,
@@ -450,12 +590,20 @@ class autoSnapDriver
                 void featureAttractionUsingFeatureEdges
                 (
                     const label iter,
-                    const bool avoidSnapProblems,
-                    const scalar featureCos,
                     const bool multiRegionFeatureSnap,
+
+                    const bool detectBaffles,
+                    const bool baffleFeaturePoints,
+                    const bool releasePoints,
+                    const bool stringFeatures,
+                    const bool avoidDiagonal,
+
+                    const scalar featureCos,
+
                     const indirectPrimitivePatch& pp,
                     const scalarField& snapDist,
                     const vectorField& nearestDisp,
+                    const vectorField& nearestNormal,
 
                     const List<List<point> >& pointFaceSurfNormals,
                     const List<List<point> >& pointFaceDisp,
@@ -465,12 +613,14 @@ class autoSnapDriver
                     vectorField& patchAttraction,
                     List<pointConstraint>& patchConstraints
                 ) const;
+
                 void preventFaceSqueeze
                 (
                     const label iter,
                     const scalar featureCos,
                     const indirectPrimitivePatch& pp,
                     const scalarField& snapDist,
+                    const vectorField& nearestAttraction,
 
                     vectorField& patchAttraction,
                     List<pointConstraint>& patchConstraints
@@ -483,15 +633,19 @@ class autoSnapDriver
                 vectorField calcNearestSurfaceFeature
                 (
                     const snapParameters& snapParams,
-                    const bool avoidSnapProblems,
+                    const bool alignMeshEdges,
                     const label iter,
                     const scalar featureCos,
                     const scalar featureAttract,
                     const scalarField& snapDist,
                     const vectorField& nearestDisp,
+                    const vectorField& nearestNormal,
                     motionSmoother& meshMover,
                     vectorField& patchAttraction,
-                    List<pointConstraint>& patchConstraints
+                    List<pointConstraint>& patchConstraints,
+
+                    DynamicList<label>& splitFaces,
+                    DynamicList<labelPair>& splits
                 ) const;
 
 
@@ -545,14 +699,6 @@ public:
                 motionSmoother&
             );
 
-            //- Get points both on patch and facezone.
-            static labelList getZoneSurfacePoints
-            (
-                const fvMesh& mesh,
-                const indirectPrimitivePatch&,
-                const word& zoneName
-            );
-
             //- Helper: calculate average cell centre per point
             static tmp<pointField> avgCellCentres
             (
@@ -575,7 +721,10 @@ public:
             //  displacement of patch points.
             static vectorField calcNearestSurface
             (
+                const bool strictRegionSnap,
                 const meshRefinement& meshRefiner,
+                const labelList& globalToMasterPatch,
+                const labelList& globalToSlavePatch,
                 const scalarField& snapDist,
                 const indirectPrimitivePatch&,
                 pointField& nearestPoint,
@@ -622,6 +771,7 @@ public:
             (
                 const dictionary& snapDict,
                 const dictionary& motionDict,
+                const bool mergePatchFaces,
                 const scalar featureCos,
                 const scalar planarAngle,
                 const snapParameters& snapParams
@@ -636,6 +786,12 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef NoRepository
+#   include "autoSnapDriverTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverTemplates.C

@@ -0,0 +1,65 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "autoSnapDriver.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class FaceList>
+Foam::labelList Foam::autoSnapDriver::getFacePoints
+(
+    const indirectPrimitivePatch& pp,
+    const FaceList& faces
+)
+{
+    // Could use PrimitivePatch & localFaces to extract points but might just
+    // as well do it ourselves.
+
+    boolList pointOnZone(pp.nPoints(), false);
+
+    forAll(faces, i)
+    {
+        const face& f = faces[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);
+}
+
+
+// ************************************************************************* //

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -124,6 +124,14 @@ Foam::layerParameters::layerParameters
         readScalar(dict.lookup("minThickness"))
     ),
     featureAngle_(readScalar(dict.lookup("featureAngle"))),
+    mergePatchFacesAngle_
+    (
+        dict.lookupOrDefault<scalar>
+        (
+            "mergePatchFacesAngle",
+            featureAngle_
+        )
+    ),
     concaveAngle_
     (
         dict.lookupOrDefault("concaveAngle", defaultConcaveAngle)

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -114,6 +114,8 @@ private:
 
         scalar featureAngle_;
 
+        scalar mergePatchFacesAngle_;
+
         scalar concaveAngle_;
 
         label nGrow_;
@@ -247,6 +249,11 @@ public:
                 return featureAngle_;
             }
 
+            scalar mergePatchFacesAngle() const
+            {
+                return mergePatchFacesAngle_;
+            }
+
             scalar concaveAngle() const
             {
                 return concaveAngle_;

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,6 +27,8 @@ License
 #include "unitConversion.H"
 #include "polyMesh.H"
 #include "globalIndex.H"
+#include "Tuple2.H"
+#include "wallPolyPatch.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -44,7 +46,15 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
         )
     ),
     nBufferLayers_(readLabel(dict.lookup("nCellsBetweenLevels"))),
-    keepPoints_(pointField(1, dict.lookup("locationInMesh"))),
+    locationsOutsideMesh_
+    (
+        dict.lookupOrDefault
+        (
+            "locationsOutsideMesh",
+            pointField(0)
+        )
+    ),
+    faceZoneControls_(dict.subOrEmptyDict("faceZoneControls")),
     allowFreeStandingZoneFaces_(dict.lookup("allowFreeStandingZoneFaces")),
     useTopologicalSnapDetection_
     (
@@ -54,8 +64,35 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
     handleSnapProblems_
     (
         dict.lookupOrDefault<Switch>("handleSnapProblems", true)
+    ),
+    interfaceRefine_
+    (
+        dict.lookupOrDefault<Switch>("interfaceRefine", true)
     )
 {
+    point locationInMesh;
+    if (dict.readIfPresent("locationInMesh", locationInMesh))
+    {
+        locationsInMesh_.append(locationInMesh);
+        zonesInMesh_.append("noneIfNotSet");// special name for no cellZone
+    }
+
+    List<Tuple2<point, word> > pointsToZone;
+    if (dict.readIfPresent("locationsInMesh", pointsToZone))
+    {
+        label nZones = locationsInMesh_.size();
+        locationsInMesh_.setSize(nZones+pointsToZone.size());
+        zonesInMesh_.setSize(locationsInMesh_.size());
+
+        forAll(pointsToZone, i)
+        {
+            locationsInMesh_[nZones] = pointsToZone[i].first();
+            zonesInMesh_[nZones] = pointsToZone[i].second();
+            nZones++;
+        }
+    }
+
+
     scalar featAngle(readScalar(dict.lookup("resolveFeatureAngle")));
 
     if (featAngle < 0 || featAngle > 180)
@@ -71,8 +108,68 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::labelList Foam::refinementParameters::findCells(const polyMesh& mesh)
-const
+Foam::dictionary Foam::refinementParameters::getZoneInfo
+(
+    const word& fzName,
+    surfaceZonesInfo::faceZoneType& faceType
+) const
+{
+    dictionary patchInfo;
+    patchInfo.add("type", wallPolyPatch::typeName);
+    faceType = surfaceZonesInfo::INTERNAL;
+
+    if (faceZoneControls_.found(fzName))
+    {
+        const dictionary& fzDict = faceZoneControls_.subDict(fzName);
+
+        if (fzDict.found("patchInfo"))
+        {
+            patchInfo = fzDict.subDict("patchInfo");
+        }
+
+        word faceTypeName;
+        if (fzDict.readIfPresent("faceType", faceTypeName))
+        {
+            faceType = surfaceZonesInfo::faceZoneTypeNames[faceTypeName];
+        }
+    }
+    return patchInfo;
+}
+
+
+Foam::labelList Foam::refinementParameters::addCellZonesToMesh
+(
+    polyMesh& mesh
+) const
+{
+    labelList zoneIDs(zonesInMesh_.size(), -1);
+    forAll(zonesInMesh_, i)
+    {
+        if
+        (
+            zonesInMesh_[i] != word::null
+         && zonesInMesh_[i] != "none"
+         && zonesInMesh_[i] != "noneIfNotSet"
+        )
+        {
+            zoneIDs[i] = surfaceZonesInfo::addCellZone
+            (
+                zonesInMesh_[i],    // name
+                labelList(0),       // addressing
+                mesh
+            );
+        }
+    }
+    return zoneIDs;
+}
+
+
+Foam::labelList Foam::refinementParameters::findCells
+(
+    const bool checkInsideMesh,
+    const polyMesh& mesh,
+    const pointField& locations
+)
 {
     // Force calculation of tet-diag decomposition (for use in findCell)
     (void)mesh.tetBasePtIs();
@@ -81,13 +178,13 @@ const
     globalIndex globalCells(mesh.nCells());
 
     // Cell label per point
-    labelList cellLabels(keepPoints_.size());
+    labelList cellLabels(locations.size());
 
-    forAll(keepPoints_, i)
+    forAll(locations, i)
     {
-        const point& keepPoint = keepPoints_[i];
+        const point& location = locations[i];
 
-        label localCellI = mesh.findCell(keepPoint);
+        label localCellI = mesh.findCell(location);
 
         label globalCellI = -1;
 
@@ -98,12 +195,13 @@ const
 
         reduce(globalCellI, maxOp<label>());
 
-        if (globalCellI == -1)
+        if (checkInsideMesh && globalCellI == -1)
         {
             FatalErrorIn
             (
-                "refinementParameters::findCells(const polyMesh&) const"
-            )   << "Point " << keepPoint
+                "refinementParameters::findCells"
+                "(const polyMesh&, const pointField&) const"
+            )   << "Point " << location
                 << " is not inside the mesh or on a face or edge." << nl
                 << "Bounding box of the mesh:" << mesh.bounds()
                 << exit(FatalError);
@@ -113,10 +211,9 @@ const
         label procI = globalCells.whichProcID(globalCellI);
         label procCellI = globalCells.toLocal(procI, globalCellI);
 
-        Info<< "Found point " << keepPoint << " in cell " << procCellI
+        Info<< "Found point " << location << " in cell " << procCellI
             << " on processor " << procI << endl;
 
-
         if (globalCells.isLocal(globalCellI))
         {
             cellLabels[i] = localCellI;
@@ -130,4 +227,48 @@ const
 }
 
 
+Foam::labelList Foam::refinementParameters::zonedLocations
+(
+    const wordList& zonesInMesh
+)
+{
+    DynamicList<label> indices(zonesInMesh.size());
+
+    forAll(zonesInMesh, i)
+    {
+        if
+        (
+            zonesInMesh[i] == word::null
+        ||  zonesInMesh[i] != "noneIfNotSet"
+        )
+        {
+            indices.append(i);
+        }
+    }
+    return indices;
+}
+
+
+Foam::labelList Foam::refinementParameters::unzonedLocations
+(
+    const wordList& zonesInMesh
+)
+{
+    DynamicList<label> indices(0);
+
+    forAll(zonesInMesh, i)
+    {
+        if
+        (
+            zonesInMesh[i] != word::null
+        &&  zonesInMesh[i] == "noneIfNotSet"
+        )
+        {
+            indices.append(i);
+        }
+    }
+    return indices;
+}
+
+
 // ************************************************************************* //

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -38,6 +38,8 @@ SourceFiles
 #include "dictionary.H"
 #include "pointField.H"
 #include "Switch.H"
+#include "wordPairHashTable.H"
+#include "surfaceZonesInfo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -73,8 +75,21 @@ class refinementParameters
         //- Number of layers between different refinement levels
         const label nBufferLayers_;
 
-        //- Areas to keep
-        const pointField keepPoints_;
+
+        // Selection of areas
+
+            //- Areas not to keep
+            const pointField locationsOutsideMesh_;
+
+            //- Areas to keep
+            pointField locationsInMesh_;
+
+            //- Region for location
+            wordList zonesInMesh_;
+
+            //- Information on how to handle faces on faceZones
+            dictionary faceZoneControls_;
+
 
         //- FaceZone faces allowed which have owner and neighbour in same
         //  cellZone?
@@ -89,6 +104,8 @@ class refinementParameters
 
         Switch handleSnapProblems_;
 
+        Switch interfaceRefine_;
+
     // Private Member Functions
 
         //- Disallow default bitwise copy construct
@@ -147,9 +164,21 @@ public:
             }
 
             //- Areas to keep
-            const pointField& keepPoints() const
+            const pointField& locationsInMesh() const
             {
-                return keepPoints_;
+                return locationsInMesh_;
+            }
+
+            //- Per area the zone name
+            const wordList& zonesInMesh() const
+            {
+                return zonesInMesh_;
+            }
+
+            //- Optional points which are checked to be outside the mesh
+            const pointField& locationsOutsideMesh() const
+            {
+                return locationsOutsideMesh_;
             }
 
             //- Are zone faces allowed only inbetween different cell zones
@@ -177,11 +206,39 @@ public:
                 return handleSnapProblems_;
             }
 
+            //- Refine cell with opposite faces with different refinement level
+            bool interfaceRefine() const
+            {
+                return interfaceRefine_;
+            }
+
 
         // Other
 
-            //- Checks that cells are in mesh. Returns cells they are in.
-            labelList findCells(const polyMesh&) const;
+            //- Get patchInfo and faceType for faceZone
+            dictionary getZoneInfo
+            (
+                const word& fzName,
+                surfaceZonesInfo::faceZoneType& faceType
+            ) const;
+
+            //- Add cellZones to mesh. Return indices of cellZones (or -1)
+            labelList addCellZonesToMesh(polyMesh&) const;
+
+            //- Checks that cells are in mesh. Returns cells (or -1) they
+            //  are in.
+            static labelList findCells
+            (
+                const bool checkInsideMesh,
+                const polyMesh&,
+                const pointField& locations
+            );
+
+            //- Extract indices of named locations (so excludes 'keepPoints')
+            static labelList zonedLocations(const wordList& zonesInMesh);
+
+            //- Extract indices of unnamed locations ('keepPoints')
+            static labelList unzonedLocations(const wordList& zonesInMesh);
 
 };
 

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/wordPairHashTable.H

@@ -0,0 +1,56 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Typedef
+    Foam::wordPairHashTable
+
+Description
+    HashTable of Pair<word>
+
+Typedef
+    Foam::wordPairHashTable
+
+Description
+    HashTable of Pair<word>
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef wordPairHashTable_H
+#define wordPairHashTable_H
+
+#include "Pair.H"
+#include "HashTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef HashTable<word, Pair<word>, FixedList<word, 2>::Hash<> >
+        wordPairHashTable;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.C

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,6 +31,7 @@ License
 Foam::snapParameters::snapParameters(const dictionary& dict)
 :
     nSmoothPatch_(readLabel(dict.lookup("nSmoothPatch"))),
+    nSmoothInternal_(dict.lookupOrDefault("nSmoothInternal", 0)),
     snapTol_(readScalar(dict.lookup("tolerance"))),
     nSmoothDispl_(readLabel(dict.lookup("nSolveIter"))),
     nSnap_(readLabel(dict.lookup("nRelaxIter"))),
@@ -44,7 +45,22 @@ Foam::snapParameters::snapParameters(const dictionary& dict)
     detectNearSurfacesSnap_
     (
         dict.lookupOrDefault("detectNearSurfacesSnap", true)
-    )
+    ),
+    strictRegionSnap_
+    (
+        dict.lookupOrDefault("strictRegionSnap", false)
+    ),
+    detectBaffles_(dict.lookupOrDefault("detectBaffles", true)),
+    baffleFeaturePoints_(dict.lookupOrDefault("baffleFeaturePoints", false)),
+    releasePoints_(dict.lookupOrDefault("releasePoints", false)),
+    stringFeatures_(dict.lookupOrDefault("stringFeatures", true)),
+    avoidDiagonal_(dict.lookupOrDefault("avoidDiagonal", false)),
+    nFaceSplitInterval_
+    (
+        dict.lookupOrDefault("nFaceSplitInterval", labelMin)
+    ),
+    concaveAngle_(dict.lookupOrDefault("concaveAngle", 45)),
+    minAreaRatio_(dict.lookupOrDefault("minAreaRatio", 0.3))
 {}
 
 

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.H

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -44,10 +44,8 @@ SourceFiles
 namespace Foam
 {
 
-// Class forward declarations
-
 /*---------------------------------------------------------------------------*\
-                           Class snapParameters Declaration
+                       Class snapParameters Declaration
 \*---------------------------------------------------------------------------*/
 
 class snapParameters
@@ -56,6 +54,8 @@ class snapParameters
 
         const label nSmoothPatch_;
 
+        const label nSmoothInternal_;
+
         const scalar snapTol_;
 
         const label nSmoothDispl_;
@@ -72,6 +72,28 @@ class snapParameters
 
         const Switch detectNearSurfacesSnap_;
 
+        const Switch strictRegionSnap_;
+
+        const Switch detectBaffles_;
+
+        const Switch baffleFeaturePoints_;
+
+        const Switch releasePoints_;
+
+        const Switch stringFeatures_;
+
+        const Switch avoidDiagonal_;
+
+
+        //- How often needs face splitting be run
+        label nFaceSplitInterval_;
+
+        //- When is angle too concave too split
+        scalar concaveAngle_;
+
+        //- When is face-split not sufficiently diagonal
+        scalar minAreaRatio_;
+
 
     // Private Member Functions
 
@@ -101,6 +123,13 @@ public:
                 return nSmoothPatch_;
             }
 
+            //- Number of internal point smoothing iterations (combined with
+            //  nSmoothPatch
+            label nSmoothInternal() const
+            {
+                return nSmoothInternal_;
+            }
+
             //- Relative distance for points to be attracted by surface
             //  feature point
             //  or edge. True distance is this factor times local
@@ -123,30 +152,90 @@ public:
                 return nSnap_;
             }
 
-            label nFeatureSnap() const
-            {
-                return nFeatureSnap_;
-            }
 
-            Switch explicitFeatureSnap() const
-            {
-                return explicitFeatureSnap_;
-            }
+            // Surface snapping specific
+
+                //- Override attraction to nearest with intersection location
+                //  at near surfaces
+                Switch detectNearSurfacesSnap() const
+                {
+                    return detectNearSurfacesSnap_;
+                }
+
+                //- Attract point to corresponding surface region only
+                Switch strictRegionSnap() const
+                {
+                    return strictRegionSnap_;
+                }
+
+
+            // Feature edge snapping specific
+
+                label nFeatureSnap() const
+                {
+                    return nFeatureSnap_;
+                }
+
+                Switch explicitFeatureSnap() const
+                {
+                    return explicitFeatureSnap_;
+                }
+
+                Switch implicitFeatureSnap() const
+                {
+                    return implicitFeatureSnap_;
+                }
+
+                Switch multiRegionFeatureSnap() const
+                {
+                    return multiRegionFeatureSnap_;
+                }
+
+                Switch detectBaffles() const
+                {
+                    return detectBaffles_;
+                }
+
+                Switch baffleFeaturePoints() const
+                {
+                    return baffleFeaturePoints_;
+                }
+
+                Switch releasePoints() const
+                {
+                    return releasePoints_;
+                }
+
+                Switch stringFeatures() const
+                {
+                    return stringFeatures_;
+                }
+
+                Switch avoidDiagonal() const
+                {
+                    return avoidDiagonal_;
+                }
+
+
+                // Face splitting
+
+                    label nFaceSplitInterval() const
+                    {
+                        return nFaceSplitInterval_;
+                    }
+
+                    scalar concaveAngle() const
+                    {
+                        return concaveAngle_;
+                    }
+
+                    scalar minAreaRatio() const
+                    {
+                        return minAreaRatio_;
+                    }
 
-            Switch implicitFeatureSnap() const
-            {
-                return implicitFeatureSnap_;
-            }
 
-            Switch multiRegionFeatureSnap() const
-            {
-                return multiRegionFeatureSnap_;
-            }
 
-            Switch detectNearSurfacesSnap() const
-            {
-                return detectNearSurfacesSnap_;
-            }
 
 };
 

src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/displacementMeshMoverMotionSolver.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -42,9 +42,6 @@ namespace Foam
 }
 
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::displacementMeshMoverMotionSolver::displacementMeshMoverMotionSolver

src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/displacementMeshMoverMotionSolver.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,8 +25,8 @@ Class
     Foam::displacementMeshMoverMotionSolver
 
 Description
-    Mesh motion solver for an fvMesh.  Based on solving the cell-centre
-    Laplacian for the motion displacement.
+    Mesh motion solver for an fvMesh. Uses externalDisplacementMeshMover
+    to do the mesh motion.
 
 SourceFiles
     displacementMeshMoverMotionSolver.C

src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/displacementMotionSolverMeshMover.C

@@ -0,0 +1,316 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "displacementMotionSolverMeshMover.H"
+#include "addToRunTimeSelectionTable.H"
+#include "pointConstraints.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(displacementMotionSolverMeshMover, 1);
+
+    addToRunTimeSelectionTable
+    (
+        externalDisplacementMeshMover,
+        displacementMotionSolverMeshMover,
+        dictionary
+    );
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+bool Foam::displacementMotionSolverMeshMover::moveMesh
+(
+    const dictionary& moveDict,
+    const label nAllowableErrors,
+    labelList& checkFaces
+)
+{
+    const label nRelaxIter = readLabel(moveDict.lookup("nRelaxIter"));
+
+    meshMover_.setDisplacementPatchFields();
+
+    Info<< typeName << " : Moving mesh ..." << endl;
+
+    scalar oldErrorReduction = -1;
+
+    bool meshOk = false;
+
+    for (label iter = 0; iter < 2*nRelaxIter; ++ iter)
+    {
+        Info<< typeName << " : Iteration " << iter << endl;
+
+        if (iter == nRelaxIter)
+        {
+            Info<< typeName
+                << " : Displacement scaling for error reduction set to 0."
+                << endl;
+            oldErrorReduction = meshMover_.setErrorReduction(0.0);
+        }
+
+        if
+        (
+            meshMover_.scaleMesh
+            (
+                checkFaces,
+                baffles_,
+                meshMover_.paramDict(),
+                moveDict,
+                true,
+                nAllowableErrors
+            )
+        )
+        {
+            Info<< typeName << " : Successfully moved mesh" << endl;
+            meshOk = true;
+            break;
+        }
+    }
+
+    if (oldErrorReduction >= 0)
+    {
+        meshMover_.setErrorReduction(oldErrorReduction);
+    }
+
+    Info<< typeName << " : Finished moving mesh ..." << endl;
+
+    return meshOk;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::displacementMotionSolverMeshMover::displacementMotionSolverMeshMover
+(
+    const dictionary& dict,
+    const List<labelPair>& baffles,
+    pointVectorField& pointDisplacement
+)
+:
+    externalDisplacementMeshMover(dict, baffles, pointDisplacement),
+
+    solverPtr_
+    (
+        displacementMotionSolver::New
+        (
+            dict.lookup("solver"),
+            pointDisplacement.mesh()(),
+            IOdictionary
+            (
+                IOobject
+                (
+                    "motionSolverDict",
+                    pointDisplacement.mesh().time().constant(),
+                    pointDisplacement.db(),
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                dict
+            ),
+            pointDisplacement,
+            pointIOField
+            (
+                IOobject
+                (
+                    "points0",
+                    pointDisplacement.mesh().time().constant(),
+                    pointDisplacement.db(),
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                pointDisplacement.mesh()().points()
+            )
+        )
+    ),
+
+    adaptPatchIDs_(getFixedValueBCs(pointDisplacement)),
+    adaptPatchPtr_(getPatch(mesh(), adaptPatchIDs_)),
+
+    scale_
+    (
+        IOobject
+        (
+            "scale",
+            pointDisplacement.time().timeName(),
+            pointDisplacement.db(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        pMesh(),
+        dimensionedScalar("scale", dimless, 1.0)
+    ),
+
+    oldPoints_(mesh().points()),
+
+    meshMover_
+    (
+        const_cast<polyMesh&>(mesh()),
+        const_cast<pointMesh&>(pMesh()),
+        adaptPatchPtr_(),
+        pointDisplacement,
+        scale_,
+        oldPoints_,
+        adaptPatchIDs_,
+        dict
+    ),
+
+    fieldSmoother_(mesh())
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::displacementMotionSolverMeshMover::~displacementMotionSolverMeshMover()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+bool Foam::displacementMotionSolverMeshMover::move
+(
+    const dictionary& moveDict,
+    const label nAllowableErrors,
+    labelList& checkFaces
+)
+{
+    // Correct and smooth the patch displacements so points next to
+    // points where the extrusion was disabled also use less extrusion.
+    // Note that this has to update the pointDisplacement boundary conditions
+    // as well, not just the internal field.
+    {
+        const label nSmoothPatchThickness = readLabel
+        (
+            moveDict.lookup("nSmoothThickness")
+        );
+
+        const word minThicknessName = word(moveDict.lookup("minThicknessName"));
+
+        scalarField zeroMinThickness;
+
+        if (minThicknessName == "none")
+        {
+            zeroMinThickness = scalarField(adaptPatchPtr_().nPoints(), 0.0);
+        }
+
+        const scalarField& minThickness =
+        (
+            (minThicknessName == "none")
+          ? zeroMinThickness
+          : mesh().lookupObject<scalarField>(minThicknessName)
+        );
+
+        const PackedBoolList isPatchMasterPoint
+        (
+            meshRefinement::getMasterPoints
+            (
+                mesh(),
+                adaptPatchPtr_().meshPoints()
+            )
+        );
+
+        const PackedBoolList isPatchMasterEdge
+        (
+            meshRefinement::getMasterEdges
+            (
+                mesh(),
+                adaptPatchPtr_().meshEdges
+                (
+                    mesh().edges(),
+                    mesh().pointEdges()
+                )
+            )
+        );
+
+        // Smooth patch displacement
+
+        vectorField displacement
+        (
+            pointDisplacement().internalField(),
+            adaptPatchPtr_().meshPoints()
+        );
+
+        fieldSmoother_.minSmoothField
+        (
+            nSmoothPatchThickness,
+            isPatchMasterPoint,
+            isPatchMasterEdge,
+            adaptPatchPtr_(),
+            minThickness,
+            displacement
+        );
+
+
+        scalar resid = 0;
+
+        forAll(displacement, patchPointI)
+        {
+            const label pointI(adaptPatchPtr_().meshPoints()[patchPointI]);
+
+            resid += mag(pointDisplacement()[pointI]-displacement[patchPointI]);
+
+            pointDisplacement()[pointI] = displacement[patchPointI];
+        }
+
+        // Take over smoothed displacements on bcs
+        meshMover_.setDisplacementPatchFields();
+    }
+
+    // Use motionSolver to calculate internal displacement
+    {
+        solverPtr_->pointDisplacement() == pointDisplacement();
+        // Force solving and constraining - just so its pointDisplacement gets
+        // the correct value
+        (void)solverPtr_->newPoints();
+        pointDisplacement() == solverPtr_->pointDisplacement();
+    }
+
+    return moveMesh(moveDict, nAllowableErrors, checkFaces);
+}
+
+
+void Foam::displacementMotionSolverMeshMover::movePoints
+(
+    const pointField& p
+)
+{
+    externalDisplacementMeshMover::movePoints(p);
+
+    // Update motion solver for new geometry
+    solverPtr_->movePoints(p);
+
+    // Update motionSmoother for new geometry (moves adaptPatchPtr_)
+    meshMover_.movePoints();
+
+    // Assume current mesh location is correct (reset oldPoints, scale)
+    meshMover_.correct();
+}
+
+
+// ************************************************************************* //

src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/displacementMotionSolverMeshMover.H

@@ -0,0 +1,158 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify i
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::displacementMotionSolverMeshMover
+
+Description
+    Quality-based under-relaxation wrapped around generic
+    displacementMotionSolver.
+
+    Example of use in layer settings in snappyHexMeshDict:
+    \verbatim
+    meshShrinker displacementMotionSolver;
+    solver displacementLaplacian;
+    displacementLaplacianCoeffs
+    {
+        diffusivity  quadratic inverseDistance 1(wall);
+    }
+    \endverbatim
+
+SourceFiles
+    displacementMotionSolverMeshMover.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef displacementMotionSolverMeshMover_H
+#define displacementMotionSolverMeshMover_H
+
+#include "externalDisplacementMeshMover.H"
+#include "displacementMotionSolver.H"
+#include "motionSmootherAlgo.H"
+#include "fieldSmoother.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+              Class displacementMotionSolverMeshMover Declaration
+\*---------------------------------------------------------------------------*/
+
+class displacementMotionSolverMeshMover
+:
+    public externalDisplacementMeshMover
+{
+    // Private Data
+
+        //- Mesh motion solver
+        autoPtr<displacementMotionSolver> solverPtr_;
+
+        //- IDs of fixedValue patches that we can modify
+        const labelList adaptPatchIDs_;
+
+        //- Combined indirect fixedValue patches that we can modify
+        autoPtr<indirectPrimitivePatch> adaptPatchPtr_;
+
+        //- Scale factor for displacemen
+        pointScalarField scale_;
+
+        //- Old point field
+        pointField oldPoints_;
+
+        //- Mesh mover algorithm
+        motionSmootherAlgo meshMover_;
+
+        //- Field smoothing
+        fieldSmoother fieldSmoother_;
+
+
+    // Private Member Functions
+
+        //- Apply the mesh mover algorithm
+        bool moveMesh
+        (
+            const dictionary& moveDict,
+            const label nAllowableErrors,
+            labelList& checkFaces
+        );
+
+
+public:
+
+    //- Runtime type information
+    TypeName("displacementMotionSolver");
+
+
+    // Constructors
+
+        //- Construct from a polyMesh and an IOdictionary
+        displacementMotionSolverMeshMover
+        (
+            const dictionary& dict,
+            const List<labelPair>& baffles,
+            pointVectorField& pointDisplacemen
+        );
+
+
+    //- Destructor
+    virtual ~displacementMotionSolverMeshMover();
+
+
+    // Member Functions
+
+        //- Move mesh using current pointDisplacement boundary values.
+        //  Return true if succesful (errors on checkFaces less than
+        //  allowable). Updates pointDisplacement.
+        virtual bool move
+        (
+            const dictionary&,
+            const label nAllowableErrors,
+            labelList& checkFaces
+        );
+
+        //- Update local data for geometry changes
+        virtual void movePoints(const pointField&);
+
+        //-  Update local data for topology changes
+        virtual void updateMesh(const mapPolyMesh&)
+        {
+            notImplemented
+            (
+                "displacementMotionSolverMeshMover::updateMesh"
+                "(const mapPolyMesh&)"
+            );
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,6 +25,7 @@ License
 
 #include "externalDisplacementMeshMover.H"
 #include "mapPolyMesh.H"
+#include "zeroFixedValuePointPatchFields.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -35,6 +36,84 @@ namespace Foam
 }
 
 
+// * * * * * * * * * * * *  Protected Member Functions * * * * * * * * * * * //
+
+Foam::labelList Foam::externalDisplacementMeshMover::getFixedValueBCs
+(
+    const pointVectorField& field
+)
+{
+    DynamicList<label> adaptPatchIDs;
+
+    forAll(field.boundaryField(), patchI)
+    {
+        const pointPatchField<vector>& patchField =
+            field.boundaryField()[patchI];
+
+        if (isA<valuePointPatchField<vector> >(patchField))
+        {
+            if (isA<zeroFixedValuePointPatchField<vector> >(patchField))
+            {
+                // Special condition of fixed boundary condition. Does not
+                // get adapted
+            }
+            else
+            {
+                adaptPatchIDs.append(patchI);
+            }
+        }
+    }
+
+    return adaptPatchIDs;
+}
+
+
+Foam::autoPtr<Foam::indirectPrimitivePatch>
+Foam::externalDisplacementMeshMover::getPatch
+(
+    const polyMesh& mesh,
+    const labelList& patchIDs
+)
+{
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    // Count faces.
+    label nFaces = 0;
+
+    forAll(patchIDs, i)
+    {
+        const polyPatch& pp = patches[patchIDs[i]];
+
+        nFaces += pp.size();
+    }
+
+    // Collect faces.
+    labelList addressing(nFaces);
+    nFaces = 0;
+
+    forAll(patchIDs, i)
+    {
+        const polyPatch& pp = patches[patchIDs[i]];
+
+        label meshFaceI = pp.start();
+
+        forAll(pp, i)
+        {
+            addressing[nFaces++] = meshFaceI++;
+        }
+    }
+
+    return autoPtr<indirectPrimitivePatch>
+    (
+        new indirectPrimitivePatch
+        (
+            IndirectList<face>(mesh.faces(), addressing),
+            mesh.points()
+        )
+    );
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::externalDisplacementMeshMover::externalDisplacementMeshMover

src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -67,6 +67,19 @@ protected:
         pointVectorField& pointDisplacement_;
 
 
+    // Protected Member functions
+
+        //- Extract fixed-value patchfields
+        static labelList getFixedValueBCs(const pointVectorField&);
+
+        //- Construct patch on selected patches
+        static autoPtr<indirectPrimitivePatch> getPatch
+        (
+            const polyMesh&,
+            const labelList&
+        );
+
+
 private:
 
     // Private Member Functions
@@ -80,6 +93,7 @@ private:
         //- Disallow default bitwise assignment
         void operator=(const externalDisplacementMeshMover&);
 
+
 public:
 
     //- Runtime type information

src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/fieldSmoother/fieldSmoother.C

@@ -0,0 +1,299 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fieldSmoother.H"
+#include "pointFields.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(fieldSmoother, 0);
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fieldSmoother::fieldSmoother(const polyMesh& mesh)
+:
+    mesh_(mesh)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::fieldSmoother::~fieldSmoother()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::fieldSmoother::smoothNormals
+(
+    const label nIter,
+    const PackedBoolList& isMeshMasterPoint,
+    const PackedBoolList& isMeshMasterEdge,
+    const labelList& fixedPoints,
+    pointVectorField& normals
+) const
+{
+    // Get smoothly varying internal normals field.
+    Info<< typeName
+        << " : Smoothing normals in interior ..." << endl;
+
+    const edgeList& edges = mesh_.edges();
+
+    // Points that do not change.
+    PackedBoolList isFixedPoint(mesh_.nPoints());
+
+    // Internal points that are fixed
+    forAll(fixedPoints, i)
+    {
+        label meshPointI = fixedPoints[i];
+        isFixedPoint.set(meshPointI, 1);
+    }
+
+    // Make sure that points that are coupled to meshPoints but not on a patch
+    // are fixed as well
+    syncTools::syncPointList(mesh_, isFixedPoint, maxEqOp<unsigned int>(), 0);
+
+
+    // Correspondence between local edges/points and mesh edges/points
+    const labelList meshPoints(identity(mesh_.nPoints()));
+
+    // Calculate inverse sum of weights
+
+    scalarField edgeWeights(mesh_.nEdges());
+    scalarField invSumWeight(meshPoints.size());
+    meshRefinement::calculateEdgeWeights
+    (
+        mesh_,
+        isMeshMasterEdge,
+        meshPoints,
+        edges,
+        edgeWeights,
+        invSumWeight
+    );
+
+    vectorField average;
+    for (label iter = 0; iter < nIter; iter++)
+    {
+        meshRefinement::weightedSum
+        (
+            mesh_,
+            isMeshMasterEdge,
+            meshPoints,
+            edges,
+            edgeWeights,
+            normals,
+            average
+        );
+        average *= invSumWeight;
+
+        // Do residual calculation every so often.
+        if ((iter % 10) == 0)
+        {
+            scalar resid = meshRefinement::gAverage
+            (
+                isMeshMasterPoint,
+                mag(normals-average)()
+            );
+            Info<< "    Iteration " << iter << "   residual " << resid << 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;
+            }
+        }
+    }
+}
+
+
+void Foam::fieldSmoother::smoothPatchNormals
+(
+    const label nIter,
+    const PackedBoolList& isPatchMasterPoint,
+    const PackedBoolList& isPatchMasterEdge,
+    const indirectPrimitivePatch& adaptPatch,
+    pointField& normals
+) const
+{
+    const edgeList& edges = adaptPatch.edges();
+    const labelList& meshPoints = adaptPatch.meshPoints();
+
+    // Get smoothly varying internal normals field.
+    Info<< typeName << " : Smoothing normals ..." << endl;
+
+    scalarField edgeWeights(edges.size());
+    scalarField invSumWeight(meshPoints.size());
+    meshRefinement::calculateEdgeWeights
+    (
+        mesh_,
+        isPatchMasterEdge,
+        meshPoints,
+        edges,
+        edgeWeights,
+        invSumWeight
+    );
+
+
+    vectorField average;
+    for (label iter = 0; iter < nIter; iter++)
+    {
+        meshRefinement::weightedSum
+        (
+            mesh_,
+            isPatchMasterEdge,
+            meshPoints,
+            edges,
+            edgeWeights,
+            normals,
+            average
+        );
+        average *= invSumWeight;
+
+        // Do residual calculation every so often.
+        if ((iter % 10) == 0)
+        {
+            scalar resid = meshRefinement::gAverage
+            (
+                isPatchMasterPoint,
+                mag(normals-average)()
+            );
+            Info<< "    Iteration " << iter << "   residual " << resid << 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;
+        }
+    }
+}
+
+
+void Foam::fieldSmoother::smoothLambdaMuDisplacement
+(
+    const label nIter,
+    const PackedBoolList& isMeshMasterPoint,
+    const PackedBoolList& isMeshMasterEdge,
+    const PackedBoolList& isToBeSmoothed,
+    vectorField& displacement
+) const
+{
+    const edgeList& edges = mesh_.edges();
+
+    // Correspondence between local edges/points and mesh edges/points
+    const labelList meshPoints(identity(mesh_.nPoints()));
+
+    // Calculate inverse sum of weights
+    scalarField edgeWeights(mesh_.nEdges());
+    scalarField invSumWeight(meshPoints.size());
+    meshRefinement::calculateEdgeWeights
+    (
+        mesh_,
+        isMeshMasterEdge,
+        meshPoints,
+        edges,
+        edgeWeights,
+        invSumWeight
+    );
+
+    // Get smoothly varying patch field.
+    Info<< typeName << " : Smoothing displacement ..." << endl;
+
+    const scalar lambda = 0.33;
+    const scalar mu = -0.34;
+
+    vectorField average;
+
+    for (label iter = 0; iter < nIter; iter++)
+    {
+        meshRefinement::weightedSum
+        (
+            mesh_,
+            isMeshMasterEdge,
+            meshPoints,
+            edges,
+            edgeWeights,
+            displacement,
+            average
+        );
+        average *= invSumWeight;
+
+        forAll(displacement, i)
+        {
+            if (isToBeSmoothed[i])
+            {
+                displacement[i] = (1-lambda)*displacement[i]+lambda*average[i];
+            }
+        }
+
+        meshRefinement::weightedSum
+        (
+            mesh_,
+            isMeshMasterEdge,
+            meshPoints,
+            edges,
+            edgeWeights,
+            displacement,
+            average
+        );
+        average *= invSumWeight;
+
+
+        forAll(displacement, i)
+        {
+            if (isToBeSmoothed[i])
+            {
+                displacement[i] = (1-mu)*displacement[i]+mu*average[i];
+            }
+        }
+
+
+        // Do residual calculation every so often.
+        if ((iter % 10) == 0)
+        {
+            scalar resid = meshRefinement::gAverage
+            (
+                isMeshMasterPoint,
+                mag(displacement-average)()
+            );
+            Info<< "    Iteration " << iter << "   residual " << resid << endl;
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/fieldSmoother/fieldSmoother.H

@@ -0,0 +1,141 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::fieldSmoother
+
+Description
+    Utility functions for mesh motion solvers
+
+SourceFiles
+    fieldSmoother.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fieldSmoother_H
+#define fieldSmoother_H
+
+#include "meshRefinement.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class fieldSmoother Declaration
+\*---------------------------------------------------------------------------*/
+
+class fieldSmoother
+{
+    // Private data
+
+        //- Reference to the poly mesh
+        const polyMesh& mesh_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        fieldSmoother(const fieldSmoother&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const fieldSmoother&);
+
+
+public:
+
+    // Run-time type information
+    TypeName("fieldSmoother");
+
+
+    // Constructors
+
+        //- Construct from a polyMesh
+        fieldSmoother(const polyMesh&);
+
+
+    //- Destructor
+    virtual ~fieldSmoother();
+
+
+    // Member Functions
+
+        //- Smooth interior normals
+        void smoothNormals
+        (
+            const label nIter,
+            const PackedBoolList& isMeshMasterPoint,
+            const PackedBoolList& isMeshMasterEdge,
+            const labelList& fixedPoints,
+            pointVectorField& normals
+        ) const;
+
+        //- Smooth patch normals
+        void smoothPatchNormals
+        (
+            const label nIter,
+            const PackedBoolList& isPatchMasterPoint,
+            const PackedBoolList& isPatchMasterEdge,
+            const indirectPrimitivePatch& adaptPatch,
+            pointField& normals
+        ) const;
+
+        //- Smooth a scalar field towards, but not beyond, a minimum value
+        template <class Type>
+        void minSmoothField
+        (
+            const label nIter,
+            const PackedBoolList& isPatchMasterPoint,
+            const PackedBoolList& isPatchMasterEdge,
+            const indirectPrimitivePatch& adaptPatch,
+            const scalarField& fieldMin,
+            Field<Type>& field
+        ) const;
+
+        //- Smooth and then un-smooth a displacement
+        void smoothLambdaMuDisplacement
+        (
+            const label nIter,
+            const PackedBoolList& isMeshMasterPoint,
+            const PackedBoolList& isMeshMasterEdge,
+            const PackedBoolList& isSmoothable,
+            vectorField& displacement
+        ) const;
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "fieldSmootherTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/fieldSmoother/fieldSmootherTemplates.C

@@ -0,0 +1,103 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template <class Type>
+void Foam::fieldSmoother::minSmoothField
+(
+    const label nIter,
+    const PackedBoolList& isPatchMasterPoint,
+    const PackedBoolList& isPatchMasterEdge,
+    const indirectPrimitivePatch& adaptPatch,
+    const scalarField& fieldMinMag,
+    Field<Type>& field
+) const
+{
+    const edgeList& edges = adaptPatch.edges();
+    const labelList& meshPoints = adaptPatch.meshPoints();
+
+    scalarField edgeWeights(edges.size());
+    scalarField invSumWeight(meshPoints.size());
+    meshRefinement::calculateEdgeWeights
+    (
+        mesh_,
+        isPatchMasterEdge,
+        meshPoints,
+        edges,
+        edgeWeights,
+        invSumWeight
+    );
+
+    // Get smoothly varying patch field.
+    Info<< typeName << " : Smoothing field ..." << endl;
+
+    for (label iter = 0; iter < nIter; iter++)
+    {
+        Field<Type> average(adaptPatch.nPoints());
+        meshRefinement::weightedSum
+        (
+            mesh_,
+            isPatchMasterEdge,
+            meshPoints,
+            edges,
+            edgeWeights,
+            field,
+            average
+        );
+        average *= invSumWeight;
+
+        // Transfer to field
+        forAll(field, pointI)
+        {
+            //full smoothing neighbours + point value
+            average[pointI] = 0.5*(field[pointI]+average[pointI]);
+
+            // perform monotonic smoothing
+            if
+            (
+                mag(average[pointI]) < mag(field[pointI])
+             && mag(average[pointI]) >= mag(fieldMinMag[pointI])
+            )
+            {
+                field[pointI] = average[pointI];
+            }
+        }
+
+        // Do residual calculation every so often.
+        if ((iter % 10) == 0)
+        {
+            scalar resid = meshRefinement::gAverage
+            (
+                isPatchMasterPoint,
+                mag(field-average)()
+            );
+            Info<< "    Iteration " << iter << "   residual " << resid << endl;
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,7 +32,7 @@ License
 #include "unitConversion.H"
 #include "PatchTools.H"
 #include "OBJstream.H"
-#include "pointData.H"
+#include "PointData.H"
 #include "zeroFixedValuePointPatchFields.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -52,241 +52,11 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-Foam::labelList Foam::medialAxisMeshMover::getFixedValueBCs
-(
-    const pointVectorField& fld
-)
-{
-    DynamicList<label> adaptPatchIDs;
-    forAll(fld.boundaryField(), patchI)
-    {
-        const pointPatchField<vector>& patchFld =
-            fld.boundaryField()[patchI];
-
-        if (isA<valuePointPatchField<vector> >(patchFld))
-        {
-            if (isA<zeroFixedValuePointPatchField<vector> >(patchFld))
-            {
-                // Special condition of fixed boundary condition. Does not
-                // get adapted
-            }
-            else
-            {
-                adaptPatchIDs.append(patchI);
-            }
-        }
-    }
-    return adaptPatchIDs;
-}
-
-
-Foam::autoPtr<Foam::indirectPrimitivePatch>
-Foam::medialAxisMeshMover::getPatch
-(
-    const polyMesh& mesh,
-    const labelList& patchIDs
-)
-{
-    const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-    // Count faces.
-    label nFaces = 0;
-
-    forAll(patchIDs, i)
-    {
-        const polyPatch& pp = patches[patchIDs[i]];
-
-        nFaces += pp.size();
-    }
-
-    // Collect faces.
-    labelList addressing(nFaces);
-    nFaces = 0;
-
-    forAll(patchIDs, i)
-    {
-        const polyPatch& pp = patches[patchIDs[i]];
-
-        label meshFaceI = pp.start();
-
-        forAll(pp, i)
-        {
-            addressing[nFaces++] = meshFaceI++;
-        }
-    }
-
-    return autoPtr<indirectPrimitivePatch>
-    (
-        new indirectPrimitivePatch
-        (
-            IndirectList<face>(mesh.faces(), addressing),
-            mesh.points()
-        )
-    );
-}
-
-
-void Foam::medialAxisMeshMover::smoothPatchNormals
-(
-    const label nSmoothDisp,
-    const PackedBoolList& isPatchMasterPoint,
-    const PackedBoolList& isPatchMasterEdge,
-    pointField& normals
-) const
-{
-    const indirectPrimitivePatch& pp = adaptPatchPtr_();
-    const edgeList& edges = pp.edges();
-    const labelList& meshPoints = pp.meshPoints();
-
-    // Get smoothly varying internal normals field.
-    Info<< typeName << " : Smoothing normals ..." << endl;
-
-    scalarField edgeWeights(edges.size());
-    scalarField invSumWeight(meshPoints.size());
-    meshRefinement::calculateEdgeWeights
-    (
-        mesh(),
-        isPatchMasterEdge,
-        meshPoints,
-        edges,
-        edgeWeights,
-        invSumWeight
-    );
-
-
-    vectorField average;
-    for (label iter = 0; iter < nSmoothDisp; iter++)
-    {
-        meshRefinement::weightedSum
-        (
-            mesh(),
-            isPatchMasterEdge,
-            meshPoints,
-            edges,
-            edgeWeights,
-            normals,
-            average
-        );
-        average *= invSumWeight;
-
-        // Do residual calculation every so often.
-        if ((iter % 10) == 0)
-        {
-            scalar resid = meshRefinement::gAverage
-            (
-                isPatchMasterPoint,
-                mag(normals-average)()
-            );
-            Info<< "    Iteration " << iter << "   residual " << resid << 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::medialAxisMeshMover::smoothNormals
-(
-    const label nSmoothDisp,
-    const PackedBoolList& isMeshMasterPoint,
-    const PackedBoolList& isMeshMasterEdge,
-    const labelList& fixedPoints,
-    pointVectorField& normals
-) const
-{
-    // Get smoothly varying internal normals field.
-    Info<< typeName
-        << " : Smoothing normals in interior ..." << endl;
-
-    const edgeList& edges = mesh().edges();
-
-    // Points that do not change.
-    PackedBoolList isFixedPoint(mesh().nPoints());
-
-    // Internal points that are fixed
-    forAll(fixedPoints, i)
-    {
-        label meshPointI = fixedPoints[i];
-        isFixedPoint.set(meshPointI, 1);
-    }
-
-    // Make sure that points that are coupled to meshPoints but not on a patch
-    // are fixed as well
-    syncTools::syncPointList(mesh(), isFixedPoint, maxEqOp<unsigned int>(), 0);
-
-
-    // Correspondence between local edges/points and mesh edges/points
-    const labelList meshPoints(identity(mesh().nPoints()));
-
-    // Calculate inverse sum of weights
-
-    scalarField edgeWeights(mesh().nEdges());
-    scalarField invSumWeight(meshPoints.size());
-    meshRefinement::calculateEdgeWeights
-    (
-        mesh(),
-        isMeshMasterEdge,
-        meshPoints,
-        edges,
-        edgeWeights,
-        invSumWeight
-    );
-
-    vectorField average;
-    for (label iter = 0; iter < nSmoothDisp; iter++)
-    {
-        meshRefinement::weightedSum
-        (
-            mesh(),
-            isMeshMasterEdge,
-            meshPoints,
-            edges,
-            edgeWeights,
-            normals,
-            average
-        );
-        average *= invSumWeight;
-
-        // Do residual calculation every so often.
-        if ((iter % 10) == 0)
-        {
-            scalar resid = meshRefinement::gAverage
-            (
-                isMeshMasterPoint,
-                mag(normals-average)()
-            );
-            Info<< "    Iteration " << iter << "   residual " << resid << 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::medialAxisMeshMover::isMaxEdge
 (
-    const List<pointData>& pointWallDist,
+    const List<PointData<vector> >& pointWallDist,
     const label edgeI,
     const scalar minCos
 ) const
@@ -331,7 +101,7 @@ bool Foam::medialAxisMeshMover::isMaxEdge
     //- Detect based on extrusion vector differing for both endpoints
     //  the idea is that e.g. a sawtooth wall can still be extruded
     //  successfully as long as it is done all to the same direction.
-    if ((pointWallDist[e[0]].v() & pointWallDist[e[1]].v()) < minCos)
+    if ((pointWallDist[e[0]].data() & pointWallDist[e[1]].data()) < minCos)
     {
         return true;
     }
@@ -439,11 +209,12 @@ void Foam::medialAxisMeshMover::update(const dictionary& coeffDict)
     pointField pointNormals(PatchTools::pointNormals(mesh(), pp));
 
     // Smooth patch normal vectors
-    smoothPatchNormals
+    fieldSmoother_.smoothPatchNormals
     (
         nSmoothSurfaceNormals,
         isPatchMasterPoint,
         isPatchMasterEdge,
+        pp,
         pointNormals
     );
 
@@ -452,7 +223,7 @@ void Foam::medialAxisMeshMover::update(const dictionary& coeffDict)
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     // Distance to wall
-    List<pointData> pointWallDist(mesh().nPoints());
+    List<PointData<vector> > pointWallDist(mesh().nPoints());
 
     // Dummy additional info for PointEdgeWave
     int dummyTrackData = 0;
@@ -461,23 +232,22 @@ void Foam::medialAxisMeshMover::update(const dictionary& coeffDict)
     // 1. Calculate distance to points where displacement is specified.
     {
         // Seed data.
-        List<pointData> wallInfo(meshPoints.size());
+        List<PointData<vector> > wallInfo(meshPoints.size());
 
         forAll(meshPoints, patchPointI)
         {
             label pointI = meshPoints[patchPointI];
-            wallInfo[patchPointI] = pointData
+            wallInfo[patchPointI] = PointData<vector>
             (
                 points[pointI],
                 0.0,
-                pointI,                       // passive scalar
                 pointNormals[patchPointI]     // surface normals
             );
         }
 
         // Do all calculations
-        List<pointData> edgeWallDist(mesh().nEdges());
-        PointEdgeWave<pointData> wallDistCalc
+        List<PointData<vector> > edgeWallDist(mesh().nEdges());
+        PointEdgeWave<PointData<vector> > wallDistCalc
         (
             mesh(),
             meshPoints,
@@ -525,11 +295,11 @@ void Foam::medialAxisMeshMover::update(const dictionary& coeffDict)
     // 2. Find points with max distance and transport information back to
     //    wall.
     {
-        List<pointData> pointMedialDist(mesh().nPoints());
-        List<pointData> edgeMedialDist(mesh().nEdges());
+        List<pointEdgePoint> pointMedialDist(mesh().nPoints());
+        List<pointEdgePoint> edgeMedialDist(mesh().nEdges());
 
         // Seed point data.
-        DynamicList<pointData> maxInfo(meshPoints.size());
+        DynamicList<pointEdgePoint> maxInfo(meshPoints.size());
         DynamicList<label> maxPoints(meshPoints.size());
 
         // 1. Medial axis points
@@ -556,12 +326,10 @@ void Foam::medialAxisMeshMover::update(const dictionary& coeffDict)
                         maxPoints.append(pointI);
                         maxInfo.append
                         (
-                            pointData
+                            pointEdgePoint
                             (
                                 points[pointI],
-                                0.0,
-                                pointI,         // passive data
-                                vector::zero    // passive data
+                                0.0
                             )
                         );
                         pointMedialDist[pointI] = maxInfo.last();
@@ -612,12 +380,10 @@ void Foam::medialAxisMeshMover::update(const dictionary& coeffDict)
                             maxPoints.append(pointI);
                             maxInfo.append
                             (
-                                pointData
+                                pointEdgePoint
                                 (
                                     medialAxisPt,   //points[pointI],
-                                    magSqr(points[pointI]-medialAxisPt),//0.0,
-                                    pointI,         // passive data
-                                    vector::zero    // passive data
+                                    magSqr(points[pointI]-medialAxisPt)//0.0
                                 )
                             );
                             pointMedialDist[pointI] = maxInfo.last();
@@ -669,12 +435,10 @@ void Foam::medialAxisMeshMover::update(const dictionary& coeffDict)
                             maxPoints.append(pointI);
                             maxInfo.append
                             (
-                                pointData
+                                pointEdgePoint
                                 (
                                     points[pointI],
-                                    0.0,
-                                    pointI,         // passive data
-                                    vector::zero    // passive data
+                                    0.0
                                 )
                             );
                             pointMedialDist[pointI] = maxInfo.last();
@@ -715,7 +479,7 @@ void Foam::medialAxisMeshMover::update(const dictionary& coeffDict)
                             // Check if angle not too large.
                             scalar cosAngle =
                             (
-                               -pointWallDist[pointI].v()
+                               -pointWallDist[pointI].data()
                               & pointNormals[i]
                             );
                             if (cosAngle > slipFeatureAngleCos)
@@ -726,12 +490,10 @@ void Foam::medialAxisMeshMover::update(const dictionary& coeffDict)
                                 maxPoints.append(pointI);
                                 maxInfo.append
                                 (
-                                    pointData
+                                    pointEdgePoint
                                     (
                                         points[pointI],
-                                        0.0,
-                                        pointI,         // passive data
-                                        vector::zero    // passive data
+                                        0.0
                                     )
                                 );
                                 pointMedialDist[pointI] = maxInfo.last();
@@ -751,7 +513,7 @@ void Foam::medialAxisMeshMover::update(const dictionary& coeffDict)
         maxPoints.shrink();
 
         // Do all calculations
-        PointEdgeWave<pointData> medialDistCalc
+        PointEdgeWave<pointEdgePoint> medialDistCalc
         (
             mesh(),
             maxPoints,
@@ -794,12 +556,12 @@ void Foam::medialAxisMeshMover::update(const dictionary& coeffDict)
         }
         else
         {
-            dispVec_[i] = pointWallDist[i].v();
+            dispVec_[i] = pointWallDist[i].data();
         }
     }
 
     // Smooth normal vectors. Do not change normals on pp.meshPoints
-    smoothNormals
+    fieldSmoother_.smoothNormals
     (
         nSmoothNormals,
         isMeshMasterPoint,
@@ -1003,79 +765,6 @@ void Foam::medialAxisMeshMover::syncPatchDisplacement
 }
 
 
-void Foam::medialAxisMeshMover::minSmoothField
-(
-    const label nSmoothDisp,
-    const PackedBoolList& isPatchMasterPoint,
-    const PackedBoolList& isPatchMasterEdge,
-    const scalarField& fieldMin,
-    scalarField& field
-) const
-{
-    const indirectPrimitivePatch& pp = adaptPatchPtr_();
-    const edgeList& edges = pp.edges();
-    const labelList& meshPoints = pp.meshPoints();
-
-    scalarField edgeWeights(edges.size());
-    scalarField invSumWeight(meshPoints.size());
-    meshRefinement::calculateEdgeWeights
-    (
-        mesh(),
-        isPatchMasterEdge,
-        meshPoints,
-        edges,
-        edgeWeights,
-        invSumWeight
-    );
-
-    // Get smoothly varying patch field.
-    Info<< typeName << " : Smoothing field ..." << endl;
-
-    for (label iter = 0; iter < nSmoothDisp; iter++)
-    {
-        scalarField average(pp.nPoints());
-        meshRefinement::weightedSum
-        (
-            mesh(),
-            isPatchMasterEdge,
-            meshPoints,
-            edges,
-            edgeWeights,
-            field,
-            average
-        );
-        average *= invSumWeight;
-
-        // 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)
-        {
-            scalar resid = meshRefinement::gAverage
-            (
-                isPatchMasterPoint,
-                mag(field-average)()
-            );
-            Info<< "    Iteration " << iter << "   residual " << resid << endl;
-        }
-    }
-}
-
 // Stop layer growth where mesh wraps around edge with a
 // large feature angle
 void Foam::medialAxisMeshMover::
@@ -1204,7 +893,8 @@ handleFeatureAngleLayerTerminations
 
     //Info<< "Added "
     //    << returnReduce(nPointCounter-nOldPointCounter, sumOp<label>())
-    //    << " point not to extrude." << endl;
+    //    << " point not to extrude due to minCos "
+    //    << minCos << endl;
 }
 
 
@@ -1226,47 +916,43 @@ void Foam::medialAxisMeshMover::findIsolatedRegions
     const labelListList& pointFaces = pp.pointFaces();
     const labelList& meshPoints = pp.meshPoints();
 
-    Info<< typeName << " : Removing isolated regions ..." << endl;
+
+    Info<< typeName << " : Removing isolated regions ..." << nl
+        << indent << "- if partially extruded faces make angle < "
+        << Foam::radToDeg(Foam::acos(minCosLayerTermination)) <<  nl;
+    if (detectExtrusionIsland)
+    {
+        Info<< indent << "- if exclusively surrounded by non-extruded points"
+            << nl;
+    }
+    else
+    {
+        Info<< indent << "- if exclusively surrounded by non-extruded faces"
+            << nl;
+    }
 
     // Keep count of number of points unextruded
     label nPointCounter = 0;
 
-
-    autoPtr<OBJstream> str;
-    if (debug)
-    {
-        str.reset
-        (
-            new OBJstream
-            (
-                mesh().time().path()
-              / "islandExcludePoints_"
-              + mesh().time().timeName()
-              + ".obj"
-            )
-        );
-        Info<< typeName
-            << " : Writing points surrounded by non-extruded points to "
-            << str().name() << endl;
-    }
-
     while (true)
     {
         // Stop layer growth where mesh wraps around edge with a
         // large feature angle
-        handleFeatureAngleLayerTerminations
-        (
-            minCosLayerTermination,
-            isPatchMasterPoint,
-            meshEdges,
+        if (minCosLayerTermination > -1)
+        {
+            handleFeatureAngleLayerTerminations
+            (
+                minCosLayerTermination,
+                isPatchMasterPoint,
+                meshEdges,
 
-            extrudeStatus,
-            patchDisp,
-            nPointCounter
-        );
-
-        syncPatchDisplacement(minThickness, patchDisp, extrudeStatus);
+                extrudeStatus,
+                patchDisp,
+                nPointCounter
+            );
 
+            syncPatchDisplacement(minThickness, patchDisp, extrudeStatus);
+        }
 
 
         // Detect either:
@@ -1389,11 +1075,6 @@ void Foam::medialAxisMeshMover::findIsolatedRegions
                 {
                     nPointCounter++;
                     nChanged++;
-
-                    if (str.valid())
-                    {
-                        str().write(pp.points()[meshPoints[patchPointI]]);
-                    }
                 }
             }
         }
@@ -1483,11 +1164,6 @@ void Foam::medialAxisMeshMover::findIsolatedRegions
                     )
                     {
                         nPointCounter++;
-
-                        if (str.valid())
-                        {
-                            str().write(pp.points()[meshPoints[f[fp]]]);
-                        }
                     }
                 }
             }
@@ -1501,98 +1177,6 @@ void Foam::medialAxisMeshMover::findIsolatedRegions
 }
 
 
-void Foam::medialAxisMeshMover::smoothLambdaMuDisplacement
-(
-    const label nSmoothDisp,
-    const PackedBoolList& isMeshMasterPoint,
-    const PackedBoolList& isMeshMasterEdge,
-    vectorField& displacement
-) const
-{
-    const edgeList& edges = mesh().edges();
-
-    // Correspondence between local edges/points and mesh edges/points
-    const labelList meshPoints(identity(mesh().nPoints()));
-
-    // Calculate inverse sum of weights
-    scalarField edgeWeights(mesh().nEdges());
-    scalarField invSumWeight(meshPoints.size());
-    meshRefinement::calculateEdgeWeights
-    (
-        mesh(),
-        isMeshMasterEdge,
-        meshPoints,
-        edges,
-        edgeWeights,
-        invSumWeight
-    );
-
-    // Get smoothly varying patch field.
-    Info<< typeName << " : Smoothing displacement ..." << endl;
-
-    const scalar lambda = 0.33;
-    const scalar mu = -0.34;
-
-    vectorField average;
-
-    for (label iter = 0; iter < nSmoothDisp; iter++)
-    {
-        meshRefinement::weightedSum
-        (
-            mesh(),
-            isMeshMasterEdge,
-            meshPoints,
-            edges,
-            edgeWeights,
-            displacement,
-            average
-        );
-        average *= invSumWeight;
-
-        forAll(displacement, i)
-        {
-            if (medialRatio_[i] > SMALL && medialRatio_[i] < 1-SMALL)
-            {
-                displacement[i] = (1-lambda)*displacement[i]+lambda*average[i];
-            }
-        }
-
-        meshRefinement::weightedSum
-        (
-            mesh(),
-            isMeshMasterEdge,
-            meshPoints,
-            edges,
-            edgeWeights,
-            displacement,
-            average
-        );
-        average *= invSumWeight;
-
-
-        forAll(displacement, i)
-        {
-            if (medialRatio_[i] > SMALL && medialRatio_[i] < 1-SMALL)
-            {
-                displacement[i] = (1-mu)*displacement[i]+mu*average[i];
-            }
-        }
-
-
-        // Do residual calculation every so often.
-        if ((iter % 10) == 0)
-        {
-            scalar resid = meshRefinement::gAverage
-            (
-                isMeshMasterPoint,
-                mag(displacement-average)()
-            );
-            Info<< "    Iteration " << iter << "   residual " << resid << endl;
-        }
-    }
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::medialAxisMeshMover::medialAxisMeshMover
@@ -1630,6 +1214,7 @@ Foam::medialAxisMeshMover::medialAxisMeshMover
         adaptPatchIDs_,
         dict
     ),
+    fieldSmoother_(mesh()),
     dispVec_
     (
         IOobject
@@ -1733,10 +1318,12 @@ void Foam::medialAxisMeshMover::calculateDisplacement
     const scalar featureAngle = readScalar(coeffDict.lookup("featureAngle"));
 
     //- Stop layer growth where mesh wraps around sharp edge
-    const scalar minCosLayerTermination = Foam::cos
+    scalar layerTerminationAngle = coeffDict.lookupOrDefault<scalar>
     (
-        degToRad(0.5*featureAngle)
+        "layerTerminationAngle",
+        0.5*featureAngle
     );
+    scalar minCosLayerTermination = Foam::cos(degToRad(layerTerminationAngle));
 
     //- Smoothing wanted patch thickness
     const label nSmoothPatchThickness = readLabel
@@ -1791,9 +1378,7 @@ void Foam::medialAxisMeshMover::calculateDisplacement
     );
 
 
-    scalarField thickness(patchDisp.size());
-
-    thickness = mag(patchDisp);
+    scalarField thickness(mag(patchDisp));
 
     forAll(thickness, patchPointI)
     {
@@ -1860,7 +1445,7 @@ void Foam::medialAxisMeshMover::calculateDisplacement
             vector n =
                 patchDisp[patchPointI]
               / (mag(patchDisp[patchPointI]) + VSMALL);
-            vector mVec = mesh().points()[pointI]-medialVec_[pointI];
+            vector mVec = medialVec_[pointI]-mesh().points()[pointI];
             mVec /= mag(mVec)+VSMALL;
             thicknessRatio *= (n&mVec);
 
@@ -1947,14 +1532,17 @@ void Foam::medialAxisMeshMover::calculateDisplacement
     }
 
 
-    // smooth layer thickness on moving patch
-    minSmoothField
+    // Smooth layer thickness on moving patch. Since some locations will have
+    // disabled the extrusion this might cause big jumps in wanted displacement
+    // for neighbouring patch points. So smooth the wanted displacement
+    // before actually trying to move the mesh.
+    fieldSmoother_.minSmoothField
     (
         nSmoothPatchThickness,
         isPatchMasterPoint,
         isPatchMasterEdge,
+        pp,
         minThickness,
-
         thickness
     );
 
@@ -1962,34 +1550,33 @@ void Foam::medialAxisMeshMover::calculateDisplacement
     // Dummy additional info for PointEdgeWave
     int dummyTrackData = 0;
 
-    List<pointData> pointWallDist(mesh().nPoints());
+    List<PointData<scalar> > 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());
+        List<PointData<scalar> > edgeWallDist(mesh().nEdges());
         labelList wallPoints(meshPoints.size());
 
         // Seed data.
-        List<pointData> wallInfo(meshPoints.size());
+        List<PointData<scalar> > wallInfo(meshPoints.size());
 
         forAll(meshPoints, patchPointI)
         {
             label pointI = meshPoints[patchPointI];
             wallPoints[patchPointI] = pointI;
-            wallInfo[patchPointI] = pointData
+            wallInfo[patchPointI] = PointData<scalar>
             (
                 points[pointI],
                 0.0,
-                thickness[patchPointI],       // transport layer thickness
-                vector::zero                  // passive vector
+                thickness[patchPointI]        // transport layer thickness
             );
         }
 
         // Do all calculations
-        PointEdgeWave<pointData> wallDistCalc
+        PointEdgeWave<PointData<scalar> > wallDistCalc
         (
             mesh(),
             wallPoints,
@@ -2016,12 +1603,12 @@ void Foam::medialAxisMeshMover::calculateDisplacement
         {
             // 1) displacement on nearest wall point, scaled by medialRatio
             //    (wall distance / medial distance)
-            // 2) pointWallDist[pointI].s() is layer thickness transported
+            // 2) pointWallDist[pointI].data() is layer thickness transported
             //    from closest wall point.
             // 3) shrink in opposite direction of addedPoints
             displacement[pointI] =
                 -medialRatio_[pointI]
-                *pointWallDist[pointI].s()
+                *pointWallDist[pointI].data()
                 *dispVec_[pointI];
         }
     }
@@ -2030,11 +1617,20 @@ void Foam::medialAxisMeshMover::calculateDisplacement
     // Smear displacement away from fixed values (medialRatio=0 or 1)
     if (nSmoothDisplacement > 0)
     {
-        smoothLambdaMuDisplacement
+        PackedBoolList isToBeSmoothed(displacement.size(), false);
+
+        forAll(displacement, i)
+        {
+            isToBeSmoothed[i] =
+                medialRatio_[i] > SMALL && medialRatio_[i] < 1-SMALL;
+        }
+
+        fieldSmoother_.smoothLambdaMuDisplacement
         (
             nSmoothDisplacement,
             isMeshMasterPoint,
             isMeshMasterEdge,
+            isToBeSmoothed,
             displacement
         );
     }
@@ -2052,8 +1648,6 @@ bool Foam::medialAxisMeshMover::shrinkMesh
     const label nSnap  = readLabel(meshQualityDict.lookup("nRelaxIter"));
 
 
-
-
     // Make sure displacement boundary conditions is uptodate with
     // internal field
     meshMover_.setDisplacementPatchFields();
@@ -2119,14 +1713,6 @@ bool Foam::medialAxisMeshMover::move
     const word minThicknessName = word(moveDict.lookup("minThicknessName"));
 
 
-    // The points have moved so before calculation update
-    // the mesh and motionSolver accordingly
-    movePoints(mesh().points());
-    //
-    //// Update any point motion bcs (e.g. timevarying)
-    //pointDisplacement_.boundaryField().updateCoeffs();
-
-
     // Extract out patch-wise displacement
     const indirectPrimitivePatch& pp = adaptPatchPtr_();
 
@@ -2180,13 +1766,10 @@ void Foam::medialAxisMeshMover::movePoints(const pointField& p)
 {
     externalDisplacementMeshMover::movePoints(p);
 
-    // Update local data for new geometry
-    adaptPatchPtr_().movePoints(p);
-
-    // Update motionSmoother for new geometry
+    // Update motionSmoother for new geometry (moves adaptPatchPtr_)
     meshMover_.movePoints();
 
-    // Assume corrent mesh location is correct
+    // Assume corrent mesh location is correct (reset oldPoints, scale)
     meshMover_.correct();
 }
 

src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2014-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,10 +28,13 @@ Description
     Mesh motion solver that uses a medial axis algorithm to work
     out a fraction between the (nearest point on a) moving surface
     and the (nearest point on a) fixed surface.
+
     This fraction is then used to scale the motion. Use
     - fixedValue on all moving patches
     - zeroFixedValue on stationary patches
     - slip on all slipping patches
+    Note that the fixedValue boundary conditions might be changed by this
+    solver to enforce consistency and a valid resulting mesh.
 
 SourceFiles
     medialAxisMeshMover.C
@@ -44,13 +47,15 @@ SourceFiles
 #include "externalDisplacementMeshMover.H"
 #include "motionSmootherAlgo.H"
 #include "autoLayerDriver.H"
+#include "fieldSmoother.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-class pointData;
+template <class DataType>
+class PointData;
 
 /*---------------------------------------------------------------------------*\
              Class medialAxisMeshMover Declaration
@@ -75,6 +80,9 @@ class medialAxisMeshMover
         //- Mesh mover algorithm
         motionSmootherAlgo meshMover_;
 
+        //- Field smoothing
+        fieldSmoother fieldSmoother_;
+
 
     // Pre-calculated medial axis information
 
@@ -95,45 +103,16 @@ class medialAxisMeshMover
 
     // Private Member Functions
 
-        //- Extract fixed-value patchfields
-        static labelList getFixedValueBCs(const pointVectorField&);
-
         //- Extract bc types. Replace fixedValue derivatives with fixedValue
         wordList getPatchFieldTypes(const pointVectorField& fld);
 
-        //- Construct patch on selected patches
-        static autoPtr<indirectPrimitivePatch> getPatch
-        (
-            const polyMesh&,
-            const labelList&
-        );
-
 
         // Calculation of medial axis information
 
-            //- Smooth normals on patch
-            void smoothPatchNormals
-            (
-                const label nSmoothDisp,
-                const PackedBoolList& isMasterPoint,
-                const PackedBoolList& isMasterEdge,
-                pointField& normals
-            ) const;
-
-            //- Smooth normals on interior
-            void smoothNormals
-            (
-                const label nSmoothDisp,
-                const PackedBoolList& isMasterPoint,
-                const PackedBoolList& isMasterEdge,
-                const labelList& fixedPoints,
-                pointVectorField& normals
-            ) const;
-
             //- Is mesh edge on a cusp of displacement
             bool isMaxEdge
             (
-                const List<pointData>& pointWallDist,
+                const List<PointData<vector> >& pointWallDist,
                 const label edgeI,
                 const scalar minCos
             ) const;
@@ -161,23 +140,6 @@ class medialAxisMeshMover
                 List<autoLayerDriver::extrudeMode>& extrudeStatus
             ) const;
 
-            void smoothLambdaMuDisplacement
-            (
-                const label nSmoothDisp,
-                const PackedBoolList& isMasterPoint,
-                const PackedBoolList& isMasterEdge,
-                vectorField& displacement
-            ) const;
-
-            void minSmoothField
-            (
-                const label nSmoothDisp,
-                const PackedBoolList& isMasterPoint,
-                const PackedBoolList& isMasterEdge,
-                const scalarField& fieldMin,
-                scalarField& field
-            ) const;
-
             //- Stop layer growth at feature edges
             void handleFeatureAngleLayerTerminations
             (

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -51,6 +51,8 @@ SourceFiles
 #include "pointFieldsFwd.H"
 #include "Tuple2.H"
 #include "pointIndexHit.H"
+#include "wordPairHashTable.H"
+#include "surfaceZonesInfo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -66,14 +68,11 @@ class refinementFeatures;
 class shellSurfaces;
 class removeCells;
 class fvMeshDistribute;
-class searchableSurface;
-class regionSplit;
-class globalIndex;
 class removePoints;
 class localPointRegion;
-
 class snapParameters;
 
+
 /*---------------------------------------------------------------------------*\
                            Class meshRefinement Declaration
 \*---------------------------------------------------------------------------*/
@@ -183,6 +182,14 @@ private:
         //- Per cc-cc vector the index of the surface hit
         labelIOList surfaceIndex_;
 
+
+        // For baffle merging
+
+            //- Original patch for baffle faces that used to be on
+            //  coupled patches
+            Map<label> faceToCoupledPatch_;
+
+
         //- User supplied face based data.
         List<Tuple2<mapType, labelList> > userFaceData_;
 
@@ -190,6 +197,15 @@ private:
         //  order changes.
         wordList meshedPatches_;
 
+        //- FaceZone to master patch name
+        HashTable<word, word> faceZoneToMasterPatch_;
+
+        //- FaceZone to slave patch name
+        HashTable<word, word> faceZoneToSlavePatch_;
+
+        //- FaceZone to method to handle faces
+        HashTable<surfaceZonesInfo::faceZoneType, word> faceZoneToType_;
+
 
     // Private Member Functions
 
@@ -216,9 +232,6 @@ private:
             pointField& neiCc
         ) const;
 
-        //- Find any intersection of surface. Store in surfaceIndex_.
-        void updateIntersections(const labelList& changedFaces);
-
         //- Remove cells. Put exposedFaces into exposedPatchIDs.
         autoPtr<mapPolyMesh> doRemoveCells
         (
@@ -369,10 +382,25 @@ private:
                 const labelList& globalToSlavePatch
             ) const;
 
+            //- Calculate intersections. Return per face -1 or the global
+            //  surface region
+            void getIntersections
+            (
+                const labelList& surfacesToTest,
+                const pointField& neiCc,
+                const labelList& testFaces,
+
+                labelList& globalRegion1,
+                labelList& globalRegion2
+            ) const;
+
             //- Determine patches for baffles
             void getBafflePatches
             (
                 const labelList& globalToMasterPatch,
+                const pointField& locationsInMesh,
+                const wordList& regionsInMesh,
+
                 const labelList& neiLevel,
                 const pointField& neiCc,
                 labelList& ownPatch,
@@ -454,7 +482,9 @@ private:
             labelList markFacesOnProblemCellsGeometric
             (
                 const snapParameters& snapParams,
-                const dictionary& motionDict
+                const dictionary& motionDict,
+                const labelList& globalToMasterPatch,
+                const labelList& globalToSlavePatch
             ) const;
 
 
@@ -494,6 +524,16 @@ private:
                 labelList& cellToZone
             ) const;
 
+            //- Finds zone per cell for cells inside region for which name
+            //  is specified.
+            void findCellZoneInsideWalk
+            (
+                const pointField& locationsInMesh,
+                const wordList& regionsInMesh,
+                const labelList& blockedFace, // per face -1 or some index >= 0
+                labelList& cellToZone
+            ) const;
+
             //- Determines cell zone from cell region information.
             bool calcRegionToZone
             (
@@ -508,7 +548,7 @@ private:
             //  marked in namedSurfaceIndex regarded as blocked.
             void findCellZoneTopo
             (
-                const point& keepPoint,
+                const pointField& locationsInMesh,
                 const labelList& namedSurfaceIndex,
                 const labelList& surfaceToCellZone,
                 labelList& cellToZone
@@ -520,6 +560,26 @@ private:
                 labelList& namedSurfaceIndex
             ) const;
 
+            //- Put cells into cellZone, faces into faceZone
+            void zonify
+            (
+                const PackedBoolList& isMasterFace,
+                const labelList& cellToZone,
+                const labelList& neiCellZone,
+                const labelList& faceToZone,
+                const boolList& meshFlipMap,
+                polyTopoChange& meshMod
+            ) const;
+
+            //- Allocate faceZoneName
+            void allocateInterRegionFaceZone
+            (
+                const label ownZone,
+                const label neiZone,
+                wordPairHashTable& zonesToFaceZone,
+                HashTable<word, labelPair, labelPair::Hash<> >&
+            ) const;
+
             //- Remove any loose standing cells
             void handleSnapProblems
             (
@@ -667,6 +727,13 @@ public:
                 return surfaceIndex_;
             }
 
+            //- For faces originating from processor faces store the original
+            //  patch
+            const Map<label>& faceToCoupledPatch() const
+            {
+                return faceToCoupledPatch_;
+            }
+
             //- Additional face data that is maintained across
             //  topo changes. Every entry is a list over all faces.
             //  Bit of a hack. Additional flag to say whether to maintain master
@@ -829,16 +896,29 @@ public:
                 const bool useTopologicalSnapDetection,
                 const bool removeEdgeConnectedCells,
                 const scalarField& perpendicularAngle,
-
-                // How to handle free-standing baffles
-                const bool mergeFreeStanding,
-                const scalar freeStandingAngle,
-
                 const dictionary& motionDict,
                 Time& runTime,
                 const labelList& globalToMasterPatch,
                 const labelList& globalToSlavePatch,
-                const point& keepPoint
+                const pointField& locationsInMesh,
+                const wordList& regionsInMesh,
+                const pointField& locationsOutsideMesh
+            );
+
+            //- Merge free-standing baffles
+            void mergeFreeStandingBaffles
+            (
+                const snapParameters& snapParams,
+                const bool useTopologicalSnapDetection,
+                const bool removeEdgeConnectedCells,
+                const scalarField& perpendicularAngle,
+                const scalar planarAngle,
+                const dictionary& motionDict,
+                Time& runTime,
+                const labelList& globalToMasterPatch,
+                const labelList& globalToSlavePatch,
+                const pointField& locationsInMesh,
+                const pointField& locationsOutsideMesh
             );
 
             //- Split off (with optional buffer layers) unreachable areas
@@ -848,7 +928,10 @@ public:
                 const label nBufferLayers,
                 const labelList& globalToMasterPatch,
                 const labelList& globalToSlavePatch,
-                const point& keepPoint
+
+                const pointField& locationsInMesh,
+                const wordList& regionsInMesh,
+                const pointField& locationsOutsideMesh
             );
 
             //- Find boundary points that connect to more than one cell
@@ -859,6 +942,16 @@ public:
             //  region and split them.
             autoPtr<mapPolyMesh> dupNonManifoldPoints();
 
+            //- Find boundary points that are on faceZones of type boundary
+            //  and duplicate them
+            autoPtr<mapPolyMesh> dupNonManifoldBoundaryPoints();
+
+            //- Merge duplicate points
+            autoPtr<mapPolyMesh> mergePoints
+            (
+                const labelList& pointToDuplicate
+            );
+
             //- Create baffle for every internal face where ownPatch != -1.
             //  External faces get repatched according to ownPatch (neiPatch
             //  should be -1 for these)
@@ -868,28 +961,55 @@ public:
                 const labelList& neiPatch
             );
 
-            //- Debug helper: check faceZones are not on processor patches
-            void checkZoneFaces() const;
-
-            //- Create baffles for faces straddling zoned surfaces. Return
-            //  baffles.
-            autoPtr<mapPolyMesh> createZoneBaffles
+            //- Get zones of given type
+            labelList getZones
             (
-                const labelList& globalToMasterPatch,
-                const labelList& globalToSlavePatch,
-                List<labelPair>&
+                const List<surfaceZonesInfo::faceZoneType>& fzTypes
+            ) const;
+
+            //- Subset baffles according to zones
+            static List<labelPair> subsetBaffles
+            (
+                const polyMesh& mesh,
+                const labelList& zoneIDs,
+                const List<labelPair>& baffles
             );
 
-            //- Merge baffles. Gets pairs of faces.
-            autoPtr<mapPolyMesh> mergeBaffles(const List<labelPair>&);
+            //- Create baffles for faces on faceZones. Return created baffles
+            //  (= pairs of faces) and corresponding faceZone
+            autoPtr<mapPolyMesh> createZoneBaffles
+            (
+                const labelList& zoneIDs,
+                List<labelPair>& baffles,
+                labelList& originatingFaceZone
+            );
+
+            //- Merge baffles. Gets pairs of faces and boundary faces to move
+            //  onto (coupled) patches
+            autoPtr<mapPolyMesh> mergeBaffles
+            (
+                const List<labelPair>&,
+                const Map<label>& faceToPatch
+            );
+
+            //- Merge all baffles on faceZones
+            autoPtr<mapPolyMesh> mergeZoneBaffles
+            (
+                const bool doInternalZones,
+                const bool doBaffleZones
+            );
 
             //- Put faces/cells into zones according to surface specification.
-            //  Returns null if no zone surfaces present. Region containing
-            //  the keepPoint will not be put into a cellZone.
+            //  Returns null if no zone surfaces present. Regions containing
+            //  locationsInMesh/regionsInMesh will be put in corresponding
+            //  cellZone. keepPoints is for backwards compatibility and sets
+            //  all yet unassigned cells to be non-zoned (zone = -1)
             autoPtr<mapPolyMesh> zonify
             (
-                const point& keepPoint,
-                const bool allowFreeStandingZoneFaces
+                const bool allowFreeStandingZoneFaces,
+                const pointField& locationsInMesh,
+                const wordList& regionsInMesh,
+                wordPairHashTable& zonesToFaceZone
             );
 
 
@@ -914,9 +1034,32 @@ public:
             //- Get patchIDs for patches added in addMeshedPatch.
             labelList meshedPatches() const;
 
+            //- Add/lookup faceZone and update information. Return index of
+            //  faceZone
+            label addFaceZone
+            (
+                const word& fzName,
+                const word& masterPatch,
+                const word& slavePatch,
+                const surfaceZonesInfo::faceZoneType& fzType
+            );
+
+            //- Lookup faceZone information. Return false if no information
+            //  for faceZone
+            bool getFaceZoneInfo
+            (
+                const word& fzName,
+                label& masterPatchID,
+                label& slavePatchID,
+                surfaceZonesInfo::faceZoneType& fzType
+            ) const;
+
             //- Select coupled faces that are not collocated
             void selectSeparatedCoupledFaces(boolList&) const;
 
+            //- Find any intersection of surface. Store in surfaceIndex_.
+            void updateIntersections(const labelList& changedFaces);
+
             //- Find region point is in. Uses optional perturbation to re-test.
             static label findRegion
             (
@@ -926,19 +1069,44 @@ public:
                 const point& p
             );
 
-            //- Split mesh. Keep part containing point.
+            static void findRegions
+            (
+                const polyMesh&,
+                const vector& perturbVec,
+                const pointField& locationsInMesh,
+                const pointField& locationsOutsideMesh,
+                const label nRegions,
+                labelList& cellRegion
+            );
+
+            //- Split mesh. Keep part containing point. Return empty map if
+            //  no cells removed.
             autoPtr<mapPolyMesh> splitMeshRegions
             (
                 const labelList& globalToMasterPatch,
                 const labelList& globalToSlavePatch,
-                const point& keepPoint
+                const pointField& locationsInMesh,
+                const pointField& locationsOutsideMesh
             );
 
             //- Split faces into two
-            autoPtr<mapPolyMesh> splitFaces
+            void doSplitFaces
             (
                 const labelList& splitFaces,
-                const labelPairList& splits
+                const labelPairList& splits,
+                polyTopoChange& meshMod
+            ) const;
+
+            //- Split faces along diagonal. Maintain mesh quality. Return
+            //  total number of faces split.
+            label splitFacesUndo
+            (
+                const labelList& splitFaces,
+                const labelPairList& splits,
+                const dictionary& motionDict,
+
+                labelList& duplicateFace,
+                List<labelPair>& baffles
             );
 
             //- Update local numbering for mesh redistribution
@@ -986,6 +1154,16 @@ public:
 
             // Merging coplanar faces and edges
 
+                //- Merge coplanar faces if sets are of size mergeSize
+                //  (usually 4)
+                label mergePatchFaces
+                (
+                    const scalar minCos,
+                    const scalar concaveCos,
+                    const label mergeSize,
+                    const labelList& patchIDs
+                );
+
                 //- Merge coplanar faces. preserveFaces is != -1 for faces
                 //  to be preserved
                 label mergePatchFacesUndo

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementMerge.C

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -34,110 +34,115 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-//// Merge faces that are in-line.
-//Foam::label Foam::meshRefinement::mergePatchFaces
-//(
-//    const scalar minCos,
-//    const scalar concaveCos,
-//    const labelList& patchIDs
-//)
-//{
-//    // Patch face merging engine
-//    combineFaces faceCombiner(mesh_);
-//
-//    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-//
-//    // Pick up all candidate cells on boundary
-//    labelHashSet boundaryCells(mesh_.nFaces()-mesh_.nInternalFaces());
-//
-//    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 mergeSets
-//    (
-//        faceCombiner.getMergeSets
-//        (
-//            minCos,
-//            concaveCos,
-//            boundaryCells
-//        )
-//    );
-//
-//    label nFaceSets = returnReduce(mergeSets.size(), sumOp<label>());
-//
-//    Info<< "mergePatchFaces : Merging " << nFaceSets
-//        << " sets of faces." << endl;
-//
-//    if (nFaceSets > 0)
-//    {
-//        // Topology changes container
-//        polyTopoChange meshMod(mesh_);
-//
-//        // Merge all faces of a set into the first face of the set. Remove
-//        // unused points.
-//        faceCombiner.setRefinement(mergeSets, 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());
-//        }
-//        else
-//        {
-//            // Delete mesh volumes. No other way to do this?
-//            mesh_.clearOut();
-//        }
-//
-//
-//        // Reset the instance for if in overwrite mode
-//        mesh_.setInstance(timeName());
-//
-//        faceCombiner.updateMesh(map);
-//
-//        // Get the kept faces that need to be recalculated.
-//        // Merging two boundary faces might shift the cell centre
-//        // (unless the faces are absolutely planar)
-//        labelHashSet retestFaces(6*mergeSets.size());
-//
-//        forAll(mergeSets, setI)
-//        {
-//            label oldMasterI = mergeSets[setI][0];
-//
-//            label faceI = map().reverseFaceMap()[oldMasterI];
-//
-//            // faceI is always uncoupled boundary face
-//            const cell& cFaces = mesh_.cells()[mesh_.faceOwner()[faceI]];
-//
-//            forAll(cFaces, i)
-//            {
-//                retestFaces.insert(cFaces[i]);
-//            }
-//        }
-//        updateMesh(map, retestFaces.toc());
-//    }
-//
-//
-//    return nFaceSets;
-//}
+// Merge faces that are in-line.
+Foam::label Foam::meshRefinement::mergePatchFaces
+(
+    const scalar minCos,
+    const scalar concaveCos,
+    const label mergeSize,
+    const labelList& patchIDs
+)
+{
+    // Patch face merging engine
+    combineFaces faceCombiner(mesh_, false);
+
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+    // Pick up all candidate cells on boundary
+    labelHashSet boundaryCells(mesh_.nFaces()-mesh_.nInternalFaces());
+
+    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 mergeSets
+    (
+        faceCombiner.getMergeSets
+        (
+            minCos,
+            concaveCos,
+            boundaryCells
+        )
+    );
+
+    if (mergeSize != -1)
+    {
+        // Keep only those that are mergeSize faces
+        label compactI = 0;
+        forAll(mergeSets, setI)
+        {
+            if (mergeSets[setI].size() == mergeSize)
+            {
+                mergeSets[compactI++] = mergeSets[setI];
+            }
+        }
+        mergeSets.setSize(compactI);
+    }
+
+
+    label nFaceSets = returnReduce(mergeSets.size(), sumOp<label>());
+
+    Info<< "Merging " << nFaceSets << " sets of faces." << nl << endl;
+
+    if (nFaceSets > 0)
+    {
+        // Topology changes container
+        polyTopoChange meshMod(mesh_);
+
+        // Merge all faces of a set into the first face of the set. Remove
+        // unused points.
+        faceCombiner.setRefinement(mergeSets, 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());
+        }
+        else
+        {
+            // Delete mesh volumes. No other way to do this?
+            mesh_.clearOut();
+        }
+
+        // Reset the instance for if in overwrite mode
+        mesh_.setInstance(timeName());
+
+        faceCombiner.updateMesh(map);
+
+        // Get the kept faces that need to be recalculated.
+        // Merging two boundary faces might shift the cell centre
+        // (unless the faces are absolutely planar)
+        labelHashSet retestFaces(2*mergeSets.size());
+
+        forAll(mergeSets, setI)
+        {
+            label oldMasterI = mergeSets[setI][0];
+            retestFaces.insert(map().reverseFaceMap()[oldMasterI]);
+        }
+        updateMesh(map, growFaceCellFace(retestFaces));
+    }
+
+    return nFaceSets;
+}
+
 //
 //
 //// Remove points not used by any face or points used by only two faces where

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -1085,7 +1085,9 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
 Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
 (
     const snapParameters& snapParams,
-    const dictionary& motionDict
+    const dictionary& motionDict,
+    const labelList& globalToMasterPatch,
+    const labelList& globalToSlavePatch
 ) const
 {
     pointField oldPoints(mesh_.points());
@@ -1162,7 +1164,10 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
         (
             autoSnapDriver::calcNearestSurface
             (
+                snapParams.strictRegionSnap(),
                 *this,
+                globalToMasterPatch,
+                globalToSlavePatch,
                 snapDist,   // attraction
                 pp,
                 nearestPoint,

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -980,6 +980,7 @@ Foam::label Foam::meshRefinement::markSurfaceRefinement
     labelList surfaceMinLevel;
     surfaces_.findHigherIntersection
     (
+        shells_,
         start,
         end,
         minLevel,

src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -33,6 +33,91 @@ License
 #include "UPtrList.H"
 #include "volumeType.H"
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::labelList Foam::refinementSurfaces::findHigherLevel
+(
+    const searchableSurface& geom,
+    const shellSurfaces& shells,
+    const List<pointIndexHit>& intersectionInfo,
+    const labelList& surfaceLevel       // current level
+) const
+{
+    // See if a cached level field available
+    labelList minLevelField;
+    geom.getField(intersectionInfo, minLevelField);
+
+
+    // Detect any uncached values and do proper search
+    labelList localLevel(surfaceLevel);
+    {
+        // Check hits:
+        // 1. cached value == -1 : store for re-testing
+        // 2. cached value != -1 : use
+        // 3. uncached : use region 0 value
+
+        DynamicList<label> retestSet;
+        label nHits = 0;
+
+        forAll(intersectionInfo, i)
+        {
+            if (intersectionInfo[i].hit())
+            {
+                nHits++;
+
+                // Check if minLevelField for this surface.
+                if (minLevelField.size())
+                {
+                    if (minLevelField[i] == -1)
+                    {
+                        retestSet.append(i);
+                    }
+                    else
+                    {
+                        localLevel[i] = max(localLevel[i], minLevelField[i]);
+                    }
+                }
+                else
+                {
+                    retestSet.append(i);
+                }
+            }
+        }
+
+        label nRetest = returnReduce(retestSet.size(), sumOp<label>());
+        if (nRetest > 0)
+        {
+            reduce(nHits, sumOp<label>());
+
+            //Info<< "Retesting " << nRetest
+            //    << " out of " << nHits
+            //    << " intersections on uncached elements on geometry "
+            //    << geom.name() << endl;
+
+            pointField samples(retestSet.size());
+            forAll(retestSet, i)
+            {
+                samples[i] = intersectionInfo[retestSet[i]].hitPoint();
+            }
+            labelList shellLevel;
+            shells.findHigherLevel
+            (
+                samples,
+                UIndirectList<label>(surfaceLevel, retestSet)(),
+                shellLevel
+            );
+            forAll(retestSet, i)
+            {
+                label sampleI = retestSet[i];
+                localLevel[sampleI] = max(localLevel[sampleI], shellLevel[i]);
+            }
+        }
+    }
+
+    return localLevel;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::refinementSurfaces::refinementSurfaces
@@ -408,19 +493,19 @@ void Foam::refinementSurfaces::setMinLevelFields
     {
         const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
 
-        // Precalculation only makes sense if there are different regions
-        // (so different refinement levels possible) and there are some
+        // Cache the refinement level (max of surface level and shell level)
+        // on a per-element basis. Only makes sense if there are lots of
         // elements. Possibly should have 'enough' elements to have fine
         // enough resolution but for now just make sure we don't catch e.g.
         // searchableBox (size=6)
-        if (geom.regions().size() > 1 && geom.globalSize() > 10)
+        if (geom.globalSize() > 10)
         {
             // Representative local coordinates and bounding sphere
             pointField ctrs;
             scalarField radiusSqr;
             geom.boundingSpheres(ctrs, radiusSqr);
 
-            labelList minLevelField(ctrs.size(), -1);
+            labelList minLevelField(ctrs.size(), 0);
             {
                 // Get the element index in a roundabout way. Problem is e.g.
                 // distributed surface where local indices differ from global
@@ -447,9 +532,78 @@ void Foam::refinementSurfaces::setMinLevelFields
             labelList shellLevel;
             shells.findHigherLevel(ctrs, minLevelField, shellLevel);
 
+
+            // In case of triangulated surfaces only cache value if triangle
+            // centre and vertices are in same shell
+            if (isA<triSurface>(geom))
+            {
+                label nUncached = 0;
+
+                // Check if points differing from ctr level
+
+                const triSurface& ts = refCast<const triSurface>(geom);
+                const pointField& points = ts.points();
+
+                // Determine minimum expected level to avoid having to
+                // test lots of points
+                labelList minPointLevel(points.size(), labelMax);
+                forAll(shellLevel, triI)
+                {
+                    const labelledTri& t = ts[triI];
+                    label level = shellLevel[triI];
+                    forAll(t, tI)
+                    {
+                        minPointLevel[t[tI]] = min(minPointLevel[t[tI]], level);
+                    }
+                }
+
+
+                // See if inside any shells with higher refinement level
+                labelList pointLevel;
+                shells.findHigherLevel(points, minPointLevel, pointLevel);
+
+
+                // See if triangle centre values differ from triangle points
+                forAll(shellLevel, triI)
+                {
+                    const labelledTri& t = ts[triI];
+                    label fLevel = shellLevel[triI];
+                    if
+                    (
+                        (pointLevel[t[0]] != fLevel)
+                     || (pointLevel[t[1]] != fLevel)
+                     || (pointLevel[t[2]] != fLevel)
+                    )
+                    {
+                        //Pout<< "Detected triangle " << t.tri(ts.points())
+                        //    << " partially inside/partially outside" << endl;
+
+                        // Mark as uncached
+                        shellLevel[triI] = -1;
+                        nUncached++;
+                    }
+                }
+
+                Info<< "For geometry " << geom.name()
+                    << " detected " << returnReduce(nUncached, sumOp<label>())
+                    << " uncached triangles out of " << geom.globalSize()
+                    << endl;
+            }
+
+
+            // Combine overall level field with current shell level. Make sure
+            // to preserve -1 (from triSurfaceMeshes with triangles partly
+            // inside/outside
             forAll(minLevelField, i)
             {
-                minLevelField[i] = max(minLevelField[i], shellLevel[i]);
+                if (min(minLevelField[i], shellLevel[i]) < 0)
+                {
+                    minLevelField[i] = -1;
+                }
+                else
+                {
+                    minLevelField[i] = max(minLevelField[i], shellLevel[i]);
+                }
             }
 
             // Store minLevelField on surface
@@ -463,6 +617,8 @@ void Foam::refinementSurfaces::setMinLevelFields
 // number.
 void Foam::refinementSurfaces::findHigherIntersection
 (
+    const shellSurfaces& shells,
+
     const pointField& start,
     const pointField& end,
     const labelList& currentLevel,   // current cell refinement level
@@ -494,41 +650,48 @@ void Foam::refinementSurfaces::findHigherIntersection
         List<pointIndexHit> intersectionInfo(start.size());
         geom.findLineAny(start, end, intersectionInfo);
 
-        // See if a cached level field available
-        labelList minLevelField;
-        geom.getField(intersectionInfo, minLevelField);
-        bool haveLevelField =
-        (
-            returnReduce(minLevelField.size(), sumOp<label>())
-          > 0
-        );
 
-        if (!haveLevelField && geom.regions().size() == 1)
+        // Surface-based refinement level
+        labelList surfaceOnlyLevel(start.size(), -1);
         {
-            minLevelField = labelList
-            (
-                intersectionInfo.size(),
-                minLevel(surfI, 0)
-            );
-            haveLevelField = true;
-        }
+            // Get per intersection the region
+            labelList region;
+            geom.getRegion(intersectionInfo, region);
 
-        if (haveLevelField)
-        {
             forAll(intersectionInfo, i)
             {
-                if
-                (
-                    intersectionInfo[i].hit()
-                 && minLevelField[i] > currentLevel[i]
-                )
+                if (intersectionInfo[i].hit())
                 {
-                    surfaces[i] = surfI;    // index of surface
-                    surfaceLevel[i] = minLevelField[i];
+                    surfaceOnlyLevel[i] = minLevel(surfI, region[i]);
                 }
             }
-            return;
         }
+
+
+        // Get shell refinement level if higher
+        const labelList localLevel
+        (
+            findHigherLevel
+            (
+                geom,
+                shells,
+                intersectionInfo,
+                surfaceOnlyLevel // starting level
+            )
+        );
+
+
+        // Combine localLevel with current level
+        forAll(localLevel, i)
+        {
+            if (localLevel[i] > currentLevel[i])
+            {
+                surfaces[i] = surfI;    // index of surface
+                surfaceLevel[i] = localLevel[i];
+            }
+        }
+
+        return;
     }
 
 
@@ -546,40 +709,48 @@ void Foam::refinementSurfaces::findHigherIntersection
         // Do intersection test
         geom.findLineAny(p0, p1, intersectionInfo);
 
-        // See if a cached level field available
-        labelList minLevelField;
-        geom.getField(intersectionInfo, minLevelField);
 
-        // Copy all hits into arguments, In-place compact misses.
-        label missI = 0;
-        forAll(intersectionInfo, i)
+        // Surface-based refinement level
+        labelList surfaceOnlyLevel(intersectionInfo.size(), -1);
         {
-            // Get the minLevel for the point
-            label minLocalLevel = -1;
+            // Get per intersection the region
+            labelList region;
+            geom.getRegion(intersectionInfo, region);
 
-            if (intersectionInfo[i].hit())
+            forAll(intersectionInfo, i)
             {
-                // Check if minLevelField for this surface.
-                if (minLevelField.size())
+                if (intersectionInfo[i].hit())
                 {
-                    minLocalLevel = minLevelField[i];
-                }
-                else
-                {
-                    // Use the min level for the surface instead. Assume
-                    // single region 0.
-                    minLocalLevel = minLevel(surfI, 0);
+                    surfaceOnlyLevel[i] = minLevel(surfI, region[i]);
                 }
             }
+        }
 
 
+        // Get shell refinement level if higher
+        const labelList localLevel
+        (
+            findHigherLevel
+            (
+                geom,
+                shells,
+                intersectionInfo,
+                surfaceOnlyLevel
+            )
+        );
+
+
+        // Combine localLevel with current level
+        label missI = 0;
+        forAll(localLevel, i)
+        {
             label pointI = intersectionToPoint[i];
 
-            if (minLocalLevel > currentLevel[pointI])
+            if (localLevel[i] > currentLevel[pointI])
             {
                 // Mark point for refinement
                 surfaces[pointI] = surfI;
-                surfaceLevel[pointI] = minLocalLevel;
+                surfaceLevel[pointI] = localLevel[i];
             }
             else
             {
@@ -590,6 +761,7 @@ void Foam::refinementSurfaces::findHigherIntersection
             }
         }
 
+
         // All done? Note that this decision should be synchronised
         if (returnReduce(missI, sumOp<label>()) == 0)
         {
@@ -1061,7 +1233,7 @@ void Foam::refinementSurfaces::findNearest
 (
     const labelList& surfacesToTest,
     const pointField& samples,
-    const  scalarField& nearestDistSqr,
+    const scalarField& nearestDistSqr,
     labelList& hitSurface,
     List<pointIndexHit>& hitInfo
 ) const
@@ -1337,4 +1509,121 @@ void Foam::refinementSurfaces::findInside
 }
 
 
+void Foam::refinementSurfaces::findNearest
+(
+    const labelList& surfacesToTest,
+    const labelListList& regions,
+
+    const pointField& samples,
+    const scalarField& nearestDistSqr,
+
+    labelList& hitSurface,
+    List<pointIndexHit>& hitInfo
+) const
+{
+    labelList geometries(UIndirectList<label>(surfaces_, surfacesToTest));
+
+    // Do the tests. Note that findNearest returns index in geometries.
+    searchableSurfacesQueries::findNearest
+    (
+        allGeometry_,
+        geometries,
+        regions,
+        samples,
+        nearestDistSqr,
+        hitSurface,
+        hitInfo
+    );
+
+    // Rework the hitSurface to be surface (i.e. index into surfaces_)
+    forAll(hitSurface, pointI)
+    {
+        if (hitSurface[pointI] != -1)
+        {
+            hitSurface[pointI] = surfacesToTest[hitSurface[pointI]];
+        }
+    }
+}
+
+
+void Foam::refinementSurfaces::findNearestRegion
+(
+    const labelList& surfacesToTest,
+    const labelListList& regions,
+
+    const pointField& samples,
+    const scalarField& nearestDistSqr,
+
+    labelList& hitSurface,
+    List<pointIndexHit>& hitInfo,
+    labelList& hitRegion,
+    vectorField& hitNormal
+) const
+{
+    labelList geometries(UIndirectList<label>(surfaces_, surfacesToTest));
+
+    // Do the tests. Note that findNearest returns index in geometries.
+    searchableSurfacesQueries::findNearest
+    (
+        allGeometry_,
+        geometries,
+        regions,
+        samples,
+        nearestDistSqr,
+        hitSurface,
+        hitInfo
+    );
+
+    // Rework the hitSurface to be surface (i.e. index into surfaces_)
+    forAll(hitSurface, pointI)
+    {
+        if (hitSurface[pointI] != -1)
+        {
+            hitSurface[pointI] = surfacesToTest[hitSurface[pointI]];
+        }
+    }
+
+    // Collect the region
+    hitRegion.setSize(hitSurface.size());
+    hitRegion = -1;
+    hitNormal.setSize(hitSurface.size());
+    hitNormal = vector::zero;
+
+    forAll(surfacesToTest, i)
+    {
+        label surfI = surfacesToTest[i];
+
+        // Collect hits for surfI
+        const labelList localIndices(findIndices(hitSurface, surfI));
+
+        List<pointIndexHit> localHits
+        (
+            UIndirectList<pointIndexHit>
+            (
+                hitInfo,
+                localIndices
+            )
+        );
+
+        // Region
+        labelList localRegion;
+        allGeometry_[surfaces_[surfI]].getRegion(localHits, localRegion);
+
+        forAll(localIndices, i)
+        {
+            hitRegion[localIndices[i]] = localRegion[i];
+        }
+
+        // Normal
+        vectorField localNormal;
+        allGeometry_[surfaces_[surfI]].getNormal(localHits, localNormal);
+
+        forAll(localIndices, i)
+        {
+            hitNormal[localIndices[i]] = localNormal[i];
+        }
+    }
+}
+
+
 // ************************************************************************* //

src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -42,6 +42,7 @@ SourceFiles
 #include "vectorList.H"
 #include "pointIndexHit.H"
 #include "surfaceZonesInfo.H"
+#include "pointList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -52,8 +53,6 @@ class searchableSurfaces;
 class shellSurfaces;
 class triSurfaceMesh;
 
-typedef List<point> pointList;
-
 /*---------------------------------------------------------------------------*\
                            Class refinementSurfaces Declaration
 \*---------------------------------------------------------------------------*/
@@ -95,6 +94,16 @@ class refinementSurfaces
 
     // Private Member Functions
 
+        //- Given intersection results with geom detect local shell refinement
+        //  level (possibly cached on triangles of geom)
+        labelList findHigherLevel
+        (
+            const searchableSurface& geom,
+            const shellSurfaces& shells,
+            const List<pointIndexHit>& intersectionInfo,
+            const labelList& surfaceLevel
+        ) const;
+
         //- Disallow default bitwise copy construct
         refinementSurfaces(const refinementSurfaces&);
 
@@ -234,6 +243,8 @@ public:
             //  Return surface number and level.
             void findHigherIntersection
             (
+                const shellSurfaces& shells,
+
                 const pointField& start,
                 const pointField& end,
                 const labelList& currentLevel,  // current cell refinement level
@@ -354,6 +365,37 @@ public:
                 const pointField& pt,
                 labelList& insideSurfaces
             ) const;
+
+            // Region wise searching
+
+                //- Find nearest point on selected regions of surfaces.
+                void findNearest
+                (
+                    const labelList& surfacesToTest,
+                    const labelListList& regions,
+
+                    const pointField& samples,
+                    const scalarField& nearestDistSqr,
+
+                    labelList& hitSurface,
+                    List<pointIndexHit>& hitInfo
+                ) const;
+
+                //- Find nearest point on selected regions of surfaces.
+                void findNearestRegion
+                (
+                    const labelList& surfacesToTest,
+                    const labelListList& regions,
+
+                    const pointField& samples,
+                    const scalarField& nearestDistSqr,
+
+                    labelList& hitSurface,
+                    List<pointIndexHit>& hitInfo,
+                    labelList& hitRegion,
+                    vectorField& hitNormal
+                ) const;
+
 };
 
 

src/mesh/autoMesh/autoHexMesh/refinementSurfaces/surfaceZonesInfo.C

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -341,6 +341,37 @@ Foam::labelList Foam::surfaceZonesInfo::getInsidePointNamedSurfaces
 }
 
 
+Foam::label Foam::surfaceZonesInfo::addCellZone
+(
+    const word& name,
+    const labelList& addressing,
+    polyMesh& mesh
+)
+{
+    cellZoneMesh& cellZones = mesh.cellZones();
+
+    label zoneI = cellZones.findZoneID(name);
+
+    if (zoneI == -1)
+    {
+        zoneI = cellZones.size();
+        cellZones.setSize(zoneI+1);
+        cellZones.set
+        (
+            zoneI,
+            new cellZone
+            (
+                name,           // name
+                addressing,     // addressing
+                zoneI,          // index
+                cellZones       // cellZoneMesh
+            )
+        );
+    }
+    return zoneI;
+}
+
+
 Foam::labelList Foam::surfaceZonesInfo::addCellZonesToMesh
 (
     const PtrList<surfaceZonesInfo>& surfList,
@@ -350,8 +381,6 @@ Foam::labelList Foam::surfaceZonesInfo::addCellZonesToMesh
 {
     labelList surfaceToCellZone(surfList.size(), -1);
 
-    cellZoneMesh& cellZones = mesh.cellZones();
-
     forAll(namedSurfaces, i)
     {
         label surfI = namedSurfaces[i];
@@ -360,24 +389,12 @@ Foam::labelList Foam::surfaceZonesInfo::addCellZonesToMesh
 
         if (cellZoneName != word::null)
         {
-            label zoneI = cellZones.findZoneID(cellZoneName);
-
-            if (zoneI == -1)
-            {
-                zoneI = cellZones.size();
-                cellZones.setSize(zoneI+1);
-                cellZones.set
-                (
-                    zoneI,
-                    new cellZone
-                    (
-                        cellZoneName,   //name
-                        labelList(0),   //addressing
-                        zoneI,          //index
-                        cellZones       //cellZoneMesh
-                    )
-                );
-            }
+            label zoneI = addCellZone
+            (
+                cellZoneName,
+                labelList(0),   // addressing
+                mesh
+            );
 
             surfaceToCellZone[surfI] = zoneI;
         }
@@ -385,7 +402,7 @@ Foam::labelList Foam::surfaceZonesInfo::addCellZonesToMesh
 
     // Check they are synced
     List<wordList> allCellZones(Pstream::nProcs());
-    allCellZones[Pstream::myProcNo()] = cellZones.names();
+    allCellZones[Pstream::myProcNo()] = mesh.cellZones().names();
     Pstream::gatherList(allCellZones);
     Pstream::scatterList(allCellZones);
 
@@ -409,6 +426,40 @@ Foam::labelList Foam::surfaceZonesInfo::addCellZonesToMesh
 }
 
 
+
+Foam::label Foam::surfaceZonesInfo::addFaceZone
+(
+    const word& name,
+    const labelList& addressing,
+    const boolList& flipMap,
+    polyMesh& mesh
+)
+{
+    faceZoneMesh& faceZones = mesh.faceZones();
+
+    label zoneI = faceZones.findZoneID(name);
+
+    if (zoneI == -1)
+    {
+        zoneI = faceZones.size();
+        faceZones.setSize(zoneI+1);
+        faceZones.set
+        (
+            zoneI,
+            new faceZone
+            (
+                name,           // name
+                addressing,     // addressing
+                flipMap,        // flipMap
+                zoneI,          // index
+                faceZones       // faceZoneMesh
+            )
+        );
+    }
+    return zoneI;
+}
+
+
 Foam::labelList Foam::surfaceZonesInfo::addFaceZonesToMesh
 (
     const PtrList<surfaceZonesInfo>& surfList,
@@ -426,25 +477,13 @@ Foam::labelList Foam::surfaceZonesInfo::addFaceZonesToMesh
 
         const word& faceZoneName = surfList[surfI].faceZoneName();
 
-        label zoneI = faceZones.findZoneID(faceZoneName);
-
-        if (zoneI == -1)
-        {
-            zoneI = faceZones.size();
-            faceZones.setSize(zoneI+1);
-            faceZones.set
-            (
-                zoneI,
-                new faceZone
-                (
-                    faceZoneName,   //name
-                    labelList(0),   //addressing
-                    boolList(0),    //flipmap
-                    zoneI,          //index
-                    faceZones       //faceZoneMesh
-                )
-            );
-        }
+        label zoneI = addFaceZone
+        (
+            faceZoneName,   //name
+            labelList(0),   //addressing
+            boolList(0),    //flipmap
+            mesh
+        );
 
         surfaceToFaceZone[surfI] = zoneI;
     }

src/mesh/autoMesh/autoHexMesh/refinementSurfaces/surfaceZonesInfo.H

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -39,6 +39,7 @@ SourceFiles
 #include "word.H"
 #include "PtrList.H"
 #include "labelList.H"
+#include "boolList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -222,6 +223,13 @@ public:
                 const PtrList<surfaceZonesInfo>& surfList
             );
 
+            static label addCellZone
+            (
+                const word& name,
+                const labelList& addressing,
+                polyMesh& mesh
+            );
+
             static labelList addCellZonesToMesh
             (
                 const PtrList<surfaceZonesInfo>& surfList,
@@ -229,6 +237,14 @@ public:
                 polyMesh& mesh
             );
 
+            static label addFaceZone
+            (
+                const word& name,
+                const labelList& addressing,
+                const boolList& flipMap,
+                polyMesh& mesh
+            );
+
             static labelList addFaceZonesToMesh
             (
                 const PtrList<surfaceZonesInfo>& surfList,

tutorials/mesh/snappyHexMesh/Allrun

@@ -7,6 +7,11 @@ cd ${0%/*} || exit 1    # Run from this directory
     ./Allrun
 )
 
+(
+    cd addLayersToFaceZone || exit
+    ./Allrun
+)
+
 exit 0
 
 # These cases are links to solver test cases and are run when the Allrun