Merge branch 'master' of ssh://hunt//home/noisy3/OpenFOAM/OpenFOAM-dev

ReleaseNotes-dev

@@ -240,11 +240,14 @@
     + =snappyHexMesh=:
       + extrude across multi-processor boundaries
       + preserve faceZones during layering
+      + combining patch faces after snapping
       + *Warning*:
             - minMedianAxisAngle angle fixed w.r.t. 17x. Set to
               90 to get same behaviour as 130 in 17x.
             - nGrow did not work in 17x. Set to 0
               to get same behaviour as 130 in 17x.
+            - nSolveIter reinstated to smooth internal mesh. Set to 0
+              to get same behaviour at 17x.
             - new minTetQuality setting (see tracking). Disable for 17x
               compatibility.
 * Post-processing
@@ -309,3 +312,8 @@
 * Other
   + compilable with =clang=
     In your prefs.sh set the WM_COMPILER to Clang
+  + dlclose error: upon exit of e.g. paraFoam you can get the error
+    Inconsistency detected by ld.so: dl-close.c: 731: _dl_close: Assertion
+    `map->l_init_called' failed!
+    This seems to happen with a mix of system libraries and a thirdParty
+    compiler.

applications/solvers/basic/potentialFoam/createFields.H

@@ -47,4 +47,10 @@
 
     label pRefCell = 0;
     scalar pRefValue = 0.0;
-    setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);
+    setRefCell
+    (
+        p,
+        potentialFlow,
+        pRefCell,
+        pRefValue
+    );

applications/solvers/basic/potentialFoam/potentialFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -42,8 +42,8 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
+    #include "readControls.H"
     #include "createFields.H"
-    #include "readSIMPLEControls.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/basic/potentialFoam/readControls.H

@@ -0,0 +1,6 @@
+    const dictionary& potentialFlow =
+        mesh.solutionDict().subDict("potentialFlow");
+
+    const int nNonOrthCorr =
+        potentialFlow.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
+

applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H

@@ -64,7 +64,7 @@
         // Further limit the reciprocal time-step
         // in the vicinity of the interface
 
-        volScalarField alpha1Bar = fvc::average(alpha1);
+        volScalarField alpha1Bar(fvc::average(alpha1));
 
         rDeltaT.dimensionedInternalField() = max
         (

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,6 +37,7 @@ Description
 #include "autoLayerDriver.H"
 #include "searchableSurfaces.H"
 #include "refinementSurfaces.H"
+#include "refinementFeatures.H"
 #include "shellSurfaces.H"
 #include "decompositionMethod.H"
 #include "fvMeshDistribute.H"
@@ -250,6 +251,20 @@ int main(int argc, char *argv[])
         << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
 
 
+    // Read feature meshes
+    // ~~~~~~~~~~~~~~~~~~~
+
+    Info<< "Reading features." << endl;
+    refinementFeatures features
+    (
+        mesh,
+        refineDict.lookup("features")
+    );
+    Info<< "Read features in = "
+        << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
+
+
+
     // Refinement engine
     // ~~~~~~~~~~~~~~~~~
 
@@ -265,6 +280,7 @@ int main(int argc, char *argv[])
         mergeDist,          // tolerance used in sorting coordinates
         overwrite,          // overwrite mesh files?
         surfaces,           // for surface intersection refinement
+        features,           // for feature edges/point based refinement
         shells              // for volume (inside/outside) refinement
     );
     Info<< "Calculated surface intersections in = "
@@ -429,13 +445,19 @@ int main(int argc, char *argv[])
 
         // Snap parameters
         snapParameters snapParams(snapDict);
+        // Temporary hack to get access to resolveFeatureAngle
+        scalar curvature;
+        {
+            refinementParameters refineParams(refineDict);
+            curvature = refineParams.curvature();
+        }
 
         if (!overwrite)
         {
             const_cast<Time&>(mesh.time())++;
         }
 
-        snapDriver.doSnap(snapDict, motionDict, snapParams);
+        snapDriver.doSnap(snapDict, motionDict, curvature, snapParams);
 
         writeMesh
         (

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict

@@ -217,10 +217,10 @@ snapControls
     //  to surface
     nSmoothPatch 3;
 
-    //- Relative distance for points to be attracted by surface feature point
-    //  or edge. True distance is this factor times local
-    //  maximum edge length.
-    tolerance 4.0;
+    //- 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)
+    tolerance 2.0;
 
     //- Number of mesh displacement relaxation iterations.
     nSolveIter 30;

applications/utilities/miscellaneous/foamFormatConvert/foamFormatConvert.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -148,6 +148,11 @@ int main(int argc, char *argv[])
 #   include "setRootCase.H"
 #   include "createTime.H"
 
+    // Make sure we do not use the master-only reading since we read
+    // fields (different per processor) as dictionaries.
+    regIOobject::fileModificationChecking = regIOobject::timeStamp;
+
+
     fileName meshDir = polyMesh::meshSubDir;
     fileName regionPrefix = "";
     word regionName = polyMesh::defaultRegion;

applications/utilities/postProcessing/sampling/sample/sampleDict

@@ -150,29 +150,29 @@ surfaces
         interpolate     true;
     }
 
-    movingWall_constant
+    walls_constant
     {
         type            patch;
-        patchName       movingWall;
+        patches         ( ".*Wall.*" );
         // Optional: whether to leave as faces (=default) or triangulate
         // triangulate     false;
     }
 
-    movingWall_interpolated
+    walls_interpolated
     {
         type            patch;
-        patchName       movingWall;
+        patches         ( ".*Wall.*" );
         interpolate     true;
         // Optional: whether to leave as faces (=default) or triangulate
         // triangulate     false;
     }
 
-    movingNearWall_interpolated
+    nearWalls_interpolated
     {
         // Sample cell values off patch. Does not need to be the near-wall
         // cell, can be arbitrarily far away.
         type            patchInternalField;
-        patchName       movingWall;
+        patches         ( ".*Wall.*" );
         distance        0.0001;
         interpolate     true;
         // Optional: whether to leave as faces (=default) or triangulate

applications/utilities/preProcessing/changeDictionary/changeDictionary.C

@@ -305,6 +305,11 @@ int main(int argc, char *argv[])
         instance = args.options()["instance"];
     }
 
+    // Make sure we do not use the master-only reading since we read
+    // fields (different per processor) as dictionaries.
+    regIOobject::fileModificationChecking = regIOobject::timeStamp;
+
+
     // Get the replacement rules from a dictionary
     IOdictionary dict
     (

applications/utilities/preProcessing/foamUpgradeCyclics/foamUpgradeCyclics.C

@@ -405,6 +405,12 @@ int main(int argc, char *argv[])
 #   include "setRootCase.H"
 #   include "createTime.H"
 
+
+    // Make sure we do not use the master-only reading since we read
+    // fields (different per processor) as dictionaries.
+    regIOobject::fileModificationChecking = regIOobject::timeStamp;
+
+
     instantList timeDirs = timeSelector::select0(runTime, args);
 
     const bool isTestRun = args.optionFound("test");

applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C

@@ -40,6 +40,7 @@ Description
 #include "treeBoundBox.H"
 #include "meshTools.H"
 #include "OFstream.H"
+#include "unitConversion.H"
 
 using namespace Foam;
 
@@ -90,6 +91,55 @@ void deleteBox
 }
 
 
+// Unmark non-manifold edges if individual triangles are not features
+void unmarkBaffles
+(
+    const triSurface& surf,
+    const scalar includedAngle,
+    List<surfaceFeatures::edgeStatus>& edgeStat
+)
+{
+    scalar minCos = Foam::cos(degToRad(180.0 - includedAngle));
+
+    const labelListList& edgeFaces = surf.edgeFaces();
+
+    forAll(edgeFaces, edgeI)
+    {
+        const labelList& eFaces = edgeFaces[edgeI];
+
+        if (eFaces.size() > 2)
+        {
+            label i0 = eFaces[0];
+            //const labelledTri& f0 = surf[i0];
+            const vector& n0 = surf.faceNormals()[i0];
+
+            //Pout<< "edge:" << edgeI << " n0:" << n0 << endl;
+
+            bool same = true;
+
+            for (label i = 1; i < eFaces.size(); i++)
+            {
+                //const labelledTri& f = surf[i];
+                const vector& n = surf.faceNormals()[eFaces[i]];
+
+                //Pout<< "    mag(n&n0): " << mag(n&n0) << endl;
+
+                if (mag(n&n0) < minCos)
+                {
+                    same = false;
+                    break;
+                }
+            }
+
+            if (same)
+            {
+                edgeStat[edgeI] = surfaceFeatures::NONE;
+            }
+        }
+    }
+}
+
+
 // Main program:
 
 int main(int argc, char *argv[])
@@ -138,6 +188,11 @@ int main(int argc, char *argv[])
         "((x0 y0 z0)(x1 y1 z1))",
         "remove edges within specified bounding box"
     );
+    argList::addBoolOption
+    (
+        "writeObj",
+        "write extendedFeatureEdgeMesh obj files"
+    );
 
 #   include "setRootCase.H"
 #   include "createTime.H"
@@ -145,6 +200,7 @@ int main(int argc, char *argv[])
     Info<< "Feature line extraction is only valid on closed manifold surfaces."
         << endl;
 
+    bool writeObj = args.optionFound("writeObj");
 
     const fileName surfFileName = args[1];
     const fileName outFileName  = args[2];
@@ -153,6 +209,8 @@ int main(int argc, char *argv[])
         << "Output feature set : " << outFileName << nl
         << endl;
 
+    fileName sFeatFileName = surfFileName.lessExt().name();
+
 
     // Read
     // ~~~~
@@ -164,7 +222,6 @@ int main(int argc, char *argv[])
     Info<< endl;
 
 
-
     // Either construct features from surface&featureangle or read set.
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -187,9 +244,9 @@ int main(int argc, char *argv[])
 
         set = surfaceFeatures(surf, includedAngle);
 
-        Info<< nl << "Writing initial features" << endl;
-        set.write("initial.fSet");
-        set.writeObj("initial");
+        // Info<< nl << "Writing initial features" << endl;
+        // set.write("initial.fSet");
+        // set.writeObj("initial");
     }
     else
     {
@@ -200,7 +257,6 @@ int main(int argc, char *argv[])
             << exit(FatalError);
     }
 
-
     Info<< nl
         << "Initial feature set:" << nl
         << "    feature points : " << set.featurePoints().size() << nl
@@ -281,11 +337,12 @@ int main(int argc, char *argv[])
     surfaceFeatures newSet(surf);
     newSet.setFromStatus(edgeStat);
 
-    Info<< endl << "Writing trimmed features to " << outFileName << endl;
-    newSet.write(outFileName);
+    Info<< endl << "Writing trimmed features to "
+        << runTime.constant()/"featureEdgeMesh"/outFileName << endl;
+    newSet.write(runTime.constant()/"featureEdgeMesh"/outFileName);
 
-    Info<< endl << "Writing edge objs." << endl;
-    newSet.writeObj("final");
+    // Info<< endl << "Writing edge objs." << endl;
+    // newSet.writeObj("final");
 
     Info<< nl
         << "Final feature set:" << nl
@@ -297,20 +354,22 @@ int main(int argc, char *argv[])
         << "        internal edges : " << newSet.nInternalEdges() << nl
         << endl;
 
-    // Extracting and writing a featureEdgeMesh
+    // Extracting and writing a extendedFeatureEdgeMesh
 
     extendedFeatureEdgeMesh feMesh
     (
         newSet,
         runTime,
-        surfFileName.lessExt().name() + ".featureEdgeMesh"
+        sFeatFileName + ".extendedFeatureEdgeMesh"
     );
 
     Info<< nl << "Writing extendedFeatureEdgeMesh to " << feMesh.objectPath()
         << endl;
 
-
-    feMesh.writeObj(surfFileName.lessExt().name());
+    if (writeObj)
+    {
+        feMesh.writeObj(surfFileName.lessExt().name());
+    }
 
     feMesh.write();
 

bin/tools/CleanFunctions

@@ -95,6 +95,7 @@ cleanCase()
         > /dev/null 2>&1 \
     )
     rm -rf VTK > /dev/null 2>&1
+    rm -f 0/cellLevel 0/pointLevel
 
     if [ -e constant/polyMesh/blockMeshDict.m4 ]
     then

etc/bashrc

@@ -62,7 +62,7 @@ foamInstall=$HOME/$WM_PROJECT
 foamCompiler=system
 
 #- Compiler:
-#    WM_COMPILER = Gcc | Gcc43 | Gcc44 | Gcc45 | Clang | Icc (Intel icc)
+#    WM_COMPILER = Gcc | Gcc43 | Gcc44 | Gcc45 | Gcc46 | Clang | Icc (Intel icc)
 export WM_COMPILER=Gcc
 unset WM_COMPILER_ARCH WM_COMPILER_LIB_ARCH
 

etc/cshrc

@@ -61,7 +61,7 @@ if ( ! $?FOAM_INST_DIR ) setenv FOAM_INST_DIR $foamInstall
 set foamCompiler=system
 
 #- Compiler:
-#    WM_COMPILER = Gcc | Gcc43 | Gcc44 | Gcc45 | Clang | Icc (Intel icc)
+#    WM_COMPILER = Gcc | Gcc43 | Gcc44 | Gcc45 | Gcc46 | Clang | Icc (Intel icc)
 setenv WM_COMPILER Gcc
 setenv WM_COMPILER_ARCH # defined but empty
 unsetenv WM_COMPILER_LIB_ARCH

etc/settings.csh

@@ -207,6 +207,13 @@ case ThirdParty:
         set gmp_version=gmp-5.0.1
         set mpfr_version=mpfr-2.4.2
         breaksw
+    case Gcc46:
+    case Gcc46++0x:
+        set gcc_version=gcc-4.6.0
+        set gmp_version=gmp-5.0.1
+        set mpfr_version=mpfr-2.4.2
+        set mpc_version=mpc-0.8.1
+        breaksw
     case Gcc45:
     case Gcc45++0x:
         set gcc_version=gcc-4.5.2
@@ -381,6 +388,10 @@ case SYSTEMOPENMPI:
     setenv PLIBS "`mpicc --showme:link`"
     set libDir=`echo "$PLIBS" | sed -e 's/.*-L\([^ ]*\).*/\1/'`
 
+    # Bit of a hack: strip off 'lib' and hope this is the path to openmpi
+    # include files and libraries.
+    setenv MPI_ARCH_PATH "${libDir:h}"
+
     if ($?FOAM_VERBOSE && $?prompt) then
         echo "Using system installed MPI:"
         echo "    compile flags : $PINC"

etc/settings.sh

@@ -227,6 +227,12 @@ OpenFOAM | ThirdParty)
         gmp_version=gmp-5.0.1
         mpfr_version=mpfr-2.4.2
         ;;
+    Gcc46 | Gcc46++0x)
+        gcc_version=gcc-4.6.0
+        gmp_version=gmp-5.0.1
+        mpfr_version=mpfr-2.4.2
+        mpc_version=mpc-0.8.1
+        ;;
     Gcc45 | Gcc45++0x)
         gcc_version=gcc-4.5.2
         gmp_version=gmp-5.0.1
@@ -402,6 +408,10 @@ SYSTEMOPENMPI)
     export PLIBS="`mpicc --showme:link`"
     libDir=`echo "$PLIBS" | sed -e 's/.*-L\([^ ]*\).*/\1/'`
 
+    # Bit of a hack: strip off 'lib' and hope this is the path to openmpi
+    # include files and libraries.
+    export MPI_ARCH_PATH="${libDir%/*}"
+
     if [ "$FOAM_VERBOSE" -a "$PS1" ]
     then
         echo "Using system installed MPI:"

src/finiteVolume/cfdTools/general/porousMedia/porousZone.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -445,6 +445,13 @@ void Foam::porousZone::writeDict(Ostream& os, bool subDict) const
         dictPtr->write(os);
     }
 
+    // thermalModel
+    if (const dictionary* dictPtr = dict_.subDictPtr("thermalModel"))
+    {
+        os  << indent << "thermalModel";
+        dictPtr->write(os);
+    }
+
     os  << decrIndent << indent << token::END_BLOCK << endl;
 }
 

src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/finiteVolume/interpolation/surfaceInterpolation/schemes/pointLinear/pointLinear.C

@@ -39,8 +39,10 @@ correction
 {
     const fvMesh& mesh = this->mesh();
 
-    GeometricField<Type, pointPatchField, pointMesh> pvf =
-        volPointInterpolation::New(mesh).interpolate(vf);
+    GeometricField<Type, pointPatchField, pointMesh> pvf
+    (
+        volPointInterpolation::New(mesh).interpolate(vf)
+    );
 
     tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tsfCorr =
         linearInterpolate(vf);

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C

@@ -607,11 +607,16 @@ void Foam::InjectionModel<CloudType>::inject(TrackData& td)
                         pPtr->rho()
                     );
 
-                // Add the new parcel
-                td.cloud().addParticle(pPtr);
-
-                massAdded += pPtr->nParticle()*pPtr->mass();
-                parcelsAdded++;
+                if (pPtr->move(td, dt))
+                {
+                    td.cloud().addParticle(pPtr);
+                    massAdded += pPtr->nParticle()*pPtr->mass();
+                    parcelsAdded++;
+                }
+                else
+                {
+                    delete pPtr;
+                }
             }
         }
     }

src/mesh/autoMesh/Make/files

@@ -4,6 +4,7 @@ autoHexMeshDriver       = $(autoHexMesh)/autoHexMeshDriver
 $(autoHexMeshDriver)/autoLayerDriver.C
 $(autoHexMeshDriver)/autoLayerDriverShrink.C
 $(autoHexMeshDriver)/autoSnapDriver.C
+$(autoHexMeshDriver)/autoSnapDriverFeature.C
 $(autoHexMeshDriver)/autoRefineDriver.C
 
 $(autoHexMeshDriver)/layerParameters/layerParameters.C
@@ -16,6 +17,7 @@ $(autoHexMesh)/meshRefinement/meshRefinement.C
 $(autoHexMesh)/meshRefinement/meshRefinementMerge.C
 $(autoHexMesh)/meshRefinement/meshRefinementProblemCells.C
 $(autoHexMesh)/meshRefinement/meshRefinementRefine.C
+$(autoHexMesh)/refinementFeatures/refinementFeatures.C
 $(autoHexMesh)/refinementSurfaces/refinementSurfaces.C
 $(autoHexMesh)/shellSurfaces/shellSurfaces.C
 $(autoHexMesh)/trackedParticle/trackedParticle.C

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

@@ -59,622 +59,6 @@ defineTypeNameAndDebug(autoLayerDriver, 0);
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-Foam::label Foam::autoLayerDriver::mergePatchFacesUndo
-(
-    const scalar minCos,
-    const scalar concaveCos,
-    const dictionary& motionDict
-)
-{
-    fvMesh& mesh = meshRefiner_.mesh();
-
-    // Patch face merging engine
-    combineFaces faceCombiner(mesh, true);
-
-    // Pick up all candidate cells on boundary
-    labelHashSet boundaryCells(mesh.nFaces()-mesh.nInternalFaces());
-
-    {
-        labelList patchIDs(meshRefiner_.meshedPatches());
-
-        const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-        forAll(patchIDs, i)
-        {
-            label patchI = patchIDs[i];
-
-            const polyPatch& patch = patches[patchI];
-
-            if (!patch.coupled())
-            {
-                forAll(patch, i)
-                {
-                    boundaryCells.insert(mesh.faceOwner()[patch.start()+i]);
-                }
-            }
-        }
-    }
-
-    // Get all sets of faces that can be merged
-    labelListList allFaceSets
-    (
-        faceCombiner.getMergeSets
-        (
-            minCos,
-            concaveCos,
-            boundaryCells
-        )
-    );
-
-    label nFaceSets = returnReduce(allFaceSets.size(), sumOp<label>());
-
-    Info<< "Merging " << nFaceSets << " sets of faces." << nl << endl;
-
-    if (nFaceSets > 0)
-    {
-        if (debug)
-        {
-            faceSet allSets(mesh, "allFaceSets", allFaceSets.size());
-            forAll(allFaceSets, setI)
-            {
-                forAll(allFaceSets[setI], i)
-                {
-                    allSets.insert(allFaceSets[setI][i]);
-                }
-            }
-            Pout<< "Writing all faces to be merged to set "
-                << allSets.objectPath() << endl;
-            allSets.instance() = mesh.time().timeName();
-            allSets.write();
-        }
-
-
-        // Topology changes container
-        polyTopoChange meshMod(mesh);
-
-        // Merge all faces of a set into the first face of the set.
-        faceCombiner.setRefinement(allFaceSets, meshMod);
-
-        // Experimental: store data for all the points that have been deleted
-        meshRefiner_.storeData
-        (
-            faceCombiner.savedPointLabels(),    // points to store
-            labelList(0),                       // faces to store
-            labelList(0)                        // cells to store
-        );
-
-        // Change the mesh (no inflation)
-        autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
-
-        // Update fields
-        mesh.updateMesh(map);
-
-        // Move mesh (since morphing does not do this)
-        if (map().hasMotionPoints())
-        {
-            mesh.movePoints(map().preMotionPoints());
-        }
-        else
-        {
-            // Delete mesh volumes.
-            mesh.clearOut();
-        }
-
-        faceCombiner.updateMesh(map);
-
-        meshRefiner_.updateMesh(map, labelList(0));
-
-
-        for (label iteration = 0; iteration < 100; iteration++)
-        {
-            Info<< nl
-                << "Undo iteration " << iteration << nl
-                << "----------------" << endl;
-
-
-            // Check mesh for errors
-            // ~~~~~~~~~~~~~~~~~~~~~
-
-            faceSet errorFaces
-            (
-                mesh,
-                "errorFaces",
-                mesh.nFaces()-mesh.nInternalFaces()
-            );
-            bool hasErrors = motionSmoother::checkMesh
-            (
-                false,  // report
-                mesh,
-                motionDict,
-                errorFaces
-            );
-
-            //if (checkEdgeConnectivity)
-            //{
-            //    Info<< "Checking edge-face connectivity (duplicate faces"
-            //        << " or non-consecutive shared vertices)" << endl;
-            //
-            //    label nOldSize = errorFaces.size();
-            //
-            //    hasErrors =
-            //        mesh.checkFaceFaces
-            //        (
-            //            false,
-            //            &errorFaces
-            //        )
-            //     || hasErrors;
-            //
-            //    Info<< "Detected additional "
-            //        <<  returnReduce
-            //            (
-            //                errorFaces.size() - nOldSize,
-            //                sumOp<label>()
-            //            )
-            //        << " faces with illegal face-face connectivity"
-            //        << endl;
-            //}
-
-            if (!hasErrors)
-            {
-                break;
-            }
-
-
-            if (debug)
-            {
-                errorFaces.instance() = mesh.time().timeName();
-                Pout<< "Writing all faces in error to faceSet "
-                    << errorFaces.objectPath() << nl << endl;
-                errorFaces.write();
-            }
-
-
-            // Check any master cells for using any of the error faces
-            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-            DynamicList<label> mastersToRestore(allFaceSets.size());
-
-            forAll(allFaceSets, setI)
-            {
-                label masterFaceI = faceCombiner.masterFace()[setI];
-
-                if (masterFaceI != -1)
-                {
-                    label masterCellII = mesh.faceOwner()[masterFaceI];
-
-                    const cell& cFaces = mesh.cells()[masterCellII];
-
-                    forAll(cFaces, i)
-                    {
-                        if (errorFaces.found(cFaces[i]))
-                        {
-                            mastersToRestore.append(masterFaceI);
-                            break;
-                        }
-                    }
-                }
-            }
-            mastersToRestore.shrink();
-
-            label nRestore = returnReduce
-            (
-                mastersToRestore.size(),
-                sumOp<label>()
-            );
-
-            Info<< "Masters that need to be restored:"
-                << nRestore << endl;
-
-            if (debug)
-            {
-                faceSet restoreSet(mesh, "mastersToRestore", mastersToRestore);
-                restoreSet.instance() = mesh.time().timeName();
-                Pout<< "Writing all " << mastersToRestore.size()
-                    << " masterfaces to be restored to set "
-                    << restoreSet.objectPath() << endl;
-                restoreSet.write();
-            }
-
-
-            if (nRestore == 0)
-            {
-                break;
-            }
-
-
-            // Undo
-            // ~~~~
-
-            // Topology changes container
-            polyTopoChange meshMod(mesh);
-
-            // Merge all faces of a set into the first face of the set.
-            // Experimental:mark all points/faces/cells that have been restored.
-            Map<label> restoredPoints(0);
-            Map<label> restoredFaces(0);
-            Map<label> restoredCells(0);
-
-            faceCombiner.setUnrefinement
-            (
-                mastersToRestore,
-                meshMod,
-                restoredPoints,
-                restoredFaces,
-                restoredCells
-            );
-
-            // Change the mesh (no inflation)
-            autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
-
-            // Update fields
-            mesh.updateMesh(map);
-
-            // Move mesh (since morphing does not do this)
-            if (map().hasMotionPoints())
-            {
-                mesh.movePoints(map().preMotionPoints());
-            }
-            else
-            {
-                // Delete mesh volumes.
-                mesh.clearOut();
-            }
-
-            faceCombiner.updateMesh(map);
-
-            // Renumber restore maps
-            inplaceMapKey(map().reversePointMap(), restoredPoints);
-            inplaceMapKey(map().reverseFaceMap(), restoredFaces);
-            inplaceMapKey(map().reverseCellMap(), restoredCells);
-
-            // Experimental:restore all points/face/cells in maps
-            meshRefiner_.updateMesh
-            (
-                map,
-                labelList(0),       // changedFaces
-                restoredPoints,
-                restoredFaces,
-                restoredCells
-            );
-
-            Info<< endl;
-        }
-
-        if (debug)
-        {
-            Pout<< "Writing merged-faces mesh to time "
-                << meshRefiner_.timeName() << nl << endl;
-            mesh.write();
-        }
-    }
-    else
-    {
-        Info<< "No faces merged ..." << endl;
-    }
-
-    return nFaceSets;
-}
-
-
-// Remove points. pointCanBeDeleted is parallel synchronised.
-Foam::autoPtr<Foam::mapPolyMesh> Foam::autoLayerDriver::doRemovePoints
-(
-    removePoints& pointRemover,
-    const boolList& pointCanBeDeleted
-)
-{
-    fvMesh& mesh = meshRefiner_.mesh();
-
-    // Topology changes container
-    polyTopoChange meshMod(mesh);
-
-    pointRemover.setRefinement(pointCanBeDeleted, meshMod);
-
-    // Change the mesh (no inflation)
-    autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
-
-    // Update fields
-    mesh.updateMesh(map);
-
-    // Move mesh (since morphing does not do this)
-    if (map().hasMotionPoints())
-    {
-        mesh.movePoints(map().preMotionPoints());
-    }
-    else
-    {
-        // Delete mesh volumes.
-        mesh.clearOut();
-    }
-
-    pointRemover.updateMesh(map);
-    meshRefiner_.updateMesh(map, labelList(0));
-
-    return map;
-}
-
-
-// Restore faces (which contain removed points)
-Foam::autoPtr<Foam::mapPolyMesh> Foam::autoLayerDriver::doRestorePoints
-(
-    removePoints& pointRemover,
-    const labelList& facesToRestore
-)
-{
-    fvMesh& mesh = meshRefiner_.mesh();
-
-    // Topology changes container
-    polyTopoChange meshMod(mesh);
-
-    // Determine sets of points and faces to restore
-    labelList localFaces, localPoints;
-    pointRemover.getUnrefimentSet
-    (
-        facesToRestore,
-        localFaces,
-        localPoints
-    );
-
-    // Undo the changes on the faces that are in error.
-    pointRemover.setUnrefinement
-    (
-        localFaces,
-        localPoints,
-        meshMod
-    );
-
-    // Change the mesh (no inflation)
-    autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
-
-    // Update fields
-    mesh.updateMesh(map);
-
-    // Move mesh (since morphing does not do this)
-    if (map().hasMotionPoints())
-    {
-        mesh.movePoints(map().preMotionPoints());
-    }
-    else
-    {
-        // Delete mesh volumes.
-        mesh.clearOut();
-    }
-
-    pointRemover.updateMesh(map);
-    meshRefiner_.updateMesh(map, labelList(0));
-
-    return map;
-}
-
-
-// Collect all faces that are both in candidateFaces and in the set.
-// If coupled face also collects the coupled face.
-Foam::labelList Foam::autoLayerDriver::collectFaces
-(
-    const labelList& candidateFaces,
-    const labelHashSet& set
-) const
-{
-    const fvMesh& mesh = meshRefiner_.mesh();
-
-    // Has face been selected?
-    boolList selected(mesh.nFaces(), false);
-
-    forAll(candidateFaces, i)
-    {
-        label faceI = candidateFaces[i];
-
-        if (set.found(faceI))
-        {
-            selected[faceI] = true;
-        }
-    }
-    syncTools::syncFaceList
-    (
-        mesh,
-        selected,
-        orEqOp<bool>()      // combine operator
-    );
-
-    labelList selectedFaces(findIndices(selected, true));
-
-    return selectedFaces;
-}
-
-
-// Pick up faces of cells of faces in set.
-Foam::labelList Foam::autoLayerDriver::growFaceCellFace
-(
-    const labelHashSet& set
-) const
-{
-    const fvMesh& mesh = meshRefiner_.mesh();
-
-    boolList selected(mesh.nFaces(), false);
-
-    forAllConstIter(faceSet, set, iter)
-    {
-        label faceI = iter.key();
-
-        label own = mesh.faceOwner()[faceI];
-
-        const cell& ownFaces = mesh.cells()[own];
-        forAll(ownFaces, ownFaceI)
-        {
-            selected[ownFaces[ownFaceI]] = true;
-        }
-
-        if (mesh.isInternalFace(faceI))
-        {
-            label nbr = mesh.faceNeighbour()[faceI];
-
-            const cell& nbrFaces = mesh.cells()[nbr];
-            forAll(nbrFaces, nbrFaceI)
-            {
-                selected[nbrFaces[nbrFaceI]] = true;
-            }
-        }
-    }
-    syncTools::syncFaceList
-    (
-        mesh,
-        selected,
-        orEqOp<bool>()      // combine operator
-    );
-    return findIndices(selected, true);
-}
-
-
-// Remove points not used by any face or points used by only two faces where
-// the edges are in line
-Foam::label Foam::autoLayerDriver::mergeEdgesUndo
-(
-    const scalar minCos,
-    const dictionary& motionDict
-)
-{
-    fvMesh& mesh = meshRefiner_.mesh();
-
-    Info<< nl
-        << "Merging all points on surface that" << nl
-        << "- are used by only two boundary faces and" << nl
-        << "- make an angle with a cosine of more than " << minCos
-        << "." << nl << endl;
-
-    // Point removal analysis engine with undo
-    removePoints pointRemover(mesh, true);
-
-    // Count usage of points
-    boolList pointCanBeDeleted;
-    label nRemove = pointRemover.countPointUsage(minCos, pointCanBeDeleted);
-
-    if (nRemove > 0)
-    {
-        Info<< "Removing " << nRemove
-            << " straight edge points ..." << nl << endl;
-
-        // Remove points
-        // ~~~~~~~~~~~~~
-
-        doRemovePoints(pointRemover, pointCanBeDeleted);
-
-
-        for (label iteration = 0; iteration < 100; iteration++)
-        {
-            Info<< nl
-                << "Undo iteration " << iteration << nl
-                << "----------------" << endl;
-
-
-            // Check mesh for errors
-            // ~~~~~~~~~~~~~~~~~~~~~
-
-            faceSet errorFaces
-            (
-                mesh,
-                "errorFaces",
-                mesh.nFaces()-mesh.nInternalFaces()
-            );
-            bool hasErrors = motionSmoother::checkMesh
-            (
-                false,  // report
-                mesh,
-                motionDict,
-                errorFaces
-            );
-            //if (checkEdgeConnectivity)
-            //{
-            //    Info<< "Checking edge-face connectivity (duplicate faces"
-            //        << " or non-consecutive shared vertices)" << endl;
-            //
-            //    label nOldSize = errorFaces.size();
-            //
-            //    hasErrors =
-            //        mesh.checkFaceFaces
-            //        (
-            //            false,
-            //            &errorFaces
-            //        )
-            //     || hasErrors;
-            //
-            //    Info<< "Detected additional "
-            //        << returnReduce(errorFaces.size()-nOldSize,sumOp<label>())
-            //        << " faces with illegal face-face connectivity"
-            //        << endl;
-            //}
-
-            if (!hasErrors)
-            {
-                break;
-            }
-
-            if (debug)
-            {
-                errorFaces.instance() = mesh.time().timeName();
-                Pout<< "**Writing all faces in error to faceSet "
-                    << errorFaces.objectPath() << nl << endl;
-                errorFaces.write();
-            }
-
-            labelList masterErrorFaces
-            (
-                collectFaces
-                (
-                    pointRemover.savedFaceLabels(),
-                    errorFaces
-                )
-            );
-
-            label n = returnReduce(masterErrorFaces.size(), sumOp<label>());
-
-            Info<< "Detected " << n
-                << " error faces on boundaries that have been merged."
-                << " These will be restored to their original faces." << nl
-                << endl;
-
-            if (n == 0)
-            {
-                if (hasErrors)
-                {
-                    Info<< "Detected "
-                        << returnReduce(errorFaces.size(), sumOp<label>())
-                        << " error faces in mesh."
-                        << " Restoring neighbours of faces in error." << nl
-                        << endl;
-
-                    labelList expandedErrorFaces
-                    (
-                        growFaceCellFace
-                        (
-                            errorFaces
-                        )
-                    );
-
-                    doRestorePoints(pointRemover, expandedErrorFaces);
-                }
-
-                break;
-            }
-
-            doRestorePoints(pointRemover, masterErrorFaces);
-        }
-
-        if (debug)
-        {
-            Pout<< "Writing merged-edges mesh to time "
-                << meshRefiner_.timeName() << nl << endl;
-            mesh.write();
-        }
-    }
-    else
-    {
-        Info<< "No straight edges simplified and no points removed ..." << endl;
-    }
-
-    return nRemove;
-}
-
-
 // For debugging: Dump displacement to .obj files
 void Foam::autoLayerDriver::dumpDisplacement
 (
@@ -793,7 +177,7 @@ void Foam::autoLayerDriver::checkMeshManifold() const
             << " points where this happens to pointSet "
             << nonManifoldPoints.name() << endl;
 
-        nonManifoldPoints.instance() = mesh.time().timeName();
+        nonManifoldPoints.instance() = meshRefiner_.timeName();
         nonManifoldPoints.write();
     }
     Info<< endl;
@@ -1080,7 +464,6 @@ void Foam::autoLayerDriver::handleWarpedFaces
 }
 
 
-
 //// No extrusion on cells with multiple patch faces. There ususally is a reason
 //// why combinePatchFaces hasn't succeeded.
 //void Foam::autoLayerDriver::handleMultiplePatchFaces
@@ -1133,7 +516,7 @@ void Foam::autoLayerDriver::handleWarpedFaces
 //
 //    if (nMultiPatchCells > 0)
 //    {
-//        multiPatchCells.instance() = mesh.time().timeName();
+//        multiPatchCells.instance() = meshRefiner_.timeName();
 //        Info<< "Writing " << nMultiPatchCells
 //            << " cells with multiple (connected) patch faces to cellSet "
 //            << multiPatchCells.objectPath() << endl;
@@ -2547,9 +1930,15 @@ void Foam::autoLayerDriver::mergePatchFacesUndo
         << "      (0=straight, 180=fully concave)" << nl
         << endl;
 
-    label nChanged = mergePatchFacesUndo(minCos, concaveCos, motionDict);
+    label nChanged = meshRefiner_.mergePatchFacesUndo
+    (
+        minCos,
+        concaveCos,
+        meshRefiner_.meshedPatches(),
+        motionDict
+    );
 
-    nChanged += mergeEdgesUndo(minCos, motionDict);
+    nChanged += meshRefiner_.mergeEdgesUndo(minCos, motionDict);
 }
 
 
@@ -3191,6 +2580,7 @@ void Foam::autoLayerDriver::addLayers
                 "addedCells",
                 findIndices(flaggedCells, true)
             );
+            addedCellSet.instance() = meshRefiner_.timeName();
             Info<< "Writing "
                 << returnReduce(addedCellSet.size(), sumOp<label>())
                 << " added cells to cellSet "
@@ -3203,6 +2593,7 @@ void Foam::autoLayerDriver::addLayers
                 "layerFaces",
                 findIndices(flaggedCells, true)
             );
+            layerFacesSet.instance() = meshRefiner_.timeName();
             Info<< "Writing "
                 << returnReduce(layerFacesSet.size(), sumOp<label>())
                 << " faces inside added layer to faceSet "
@@ -3279,6 +2670,9 @@ void Foam::autoLayerDriver::addLayers
         mesh.clearOut();
     }
 
+    // Reset the instance for if in overwrite mode
+    mesh.setInstance(meshRefiner_.timeName());
+
     meshRefiner_.updateMesh(map, labelList(0));
 
 
@@ -3345,6 +2739,7 @@ void Foam::autoLayerDriver::addLayers
     // ~~~~~~~~~~
 
     cellSet addedCellSet(mesh, "addedCells", findIndices(flaggedCells, true));
+    addedCellSet.instance() = meshRefiner_.timeName();
     Info<< "Writing "
         << returnReduce(addedCellSet.size(), sumOp<label>())
         << " added cells to cellSet "
@@ -3352,6 +2747,7 @@ void Foam::autoLayerDriver::addLayers
     addedCellSet.write();
 
     faceSet layerFacesSet(mesh, "layerFaces", findIndices(flaggedFaces, true));
+    layerFacesSet.instance() = meshRefiner_.timeName();
     Info<< "Writing "
         << returnReduce(layerFacesSet.size(), sumOp<label>())
         << " faces inside added layer to faceSet "

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

@@ -105,46 +105,6 @@ class autoLayerDriver
 
     // Private Member Functions
 
-
-        // Face merging
-
-            //- Merge patch faces. Undo until no checkMesh errors.
-            label mergePatchFacesUndo
-            (
-                const scalar minCos,
-                const scalar concaveCos,
-                const dictionary&
-            );
-
-            //- Remove points.
-            autoPtr<mapPolyMesh> doRemovePoints
-            (
-                removePoints& pointRemover,
-                const boolList& pointCanBeDeleted
-            );
-
-            //- Restore faces (which contain removed points)
-            autoPtr<mapPolyMesh> doRestorePoints
-            (
-                removePoints& pointRemover,
-                const labelList& facesToRestore
-            );
-
-            //- Return candidateFaces that are also in set.
-            labelList collectFaces
-            (
-                const labelList& candidateFaces,
-                const labelHashSet& set
-            ) const;
-
-            //- Pick up faces of cells of faces in set.
-            labelList growFaceCellFace(const labelHashSet&) const;
-
-            //- Remove points not used by any face or points used by only
-            //  two faces where the edges are in line
-            label mergeEdgesUndo(const scalar minCos, const dictionary&);
-
-
         // Layers
 
             //- For debugging: Dump displacement to .obj files

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,8 +30,8 @@ License
 #include "cellSet.H"
 #include "syncTools.H"
 #include "refinementParameters.H"
-#include "featureEdgeMesh.H"
 #include "refinementSurfaces.H"
+#include "refinementFeatures.H"
 #include "shellSurfaces.H"
 #include "mapDistributePolyMesh.H"
 #include "unitConversion.H"
@@ -48,62 +48,6 @@ defineTypeNameAndDebug(autoRefineDriver, 0);
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-// Read explicit feature edges
-Foam::label Foam::autoRefineDriver::readFeatureEdges
-(
-    const PtrList<dictionary>& featDicts,
-    PtrList<featureEdgeMesh>& featureMeshes,
-    labelList& featureLevels
-) const
-{
-    Info<< "Reading external feature lines." << endl;
-
-    const fvMesh& mesh = meshRefiner_.mesh();
-
-    featureMeshes.setSize(featDicts.size());
-    featureLevels.setSize(featDicts.size());
-
-    forAll(featDicts, i)
-    {
-        const dictionary& dict = featDicts[i];
-
-        fileName featFileName(dict.lookup("file"));
-
-        featureMeshes.set
-        (
-            i,
-            new featureEdgeMesh
-            (
-                IOobject
-                (
-                    featFileName,                       // name
-                    //mesh.time().findInstance("triSurface", featFileName),
-                    //                                    // instance
-                    mesh.time().constant(),             // instance
-                    "triSurface",                       // local
-                    mesh.time(),                        // registry
-                    IOobject::MUST_READ,
-                    IOobject::NO_WRITE,
-                    false
-                )
-            )
-        );
-
-        featureMeshes[i].mergePoints(meshRefiner_.mergeDistance());
-        featureLevels[i] = readLabel(dict.lookup("level"));
-
-        Info<< "Refinement level " << featureLevels[i]
-            << " for all cells crossed by feature " << featFileName
-            << " (" << featureMeshes[i].points().size() << " points, "
-            << featureMeshes[i].edges().size() << " edges)." << endl;
-    }
-
-    Info<< "Read feature lines in = "
-        << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
-
-    return featureMeshes.size();
-}
-
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -128,23 +72,15 @@ Foam::autoRefineDriver::autoRefineDriver
 Foam::label Foam::autoRefineDriver::featureEdgeRefine
 (
     const refinementParameters& refineParams,
-    const PtrList<dictionary>& featDicts,
     const label maxIter,
     const label minRefine
 )
 {
     const fvMesh& mesh = meshRefiner_.mesh();
 
-    // Read explicit feature edges
-    PtrList<featureEdgeMesh> featureMeshes;
-    // Per feature the refinement level
-    labelList featureLevels;
-    readFeatureEdges(featDicts, featureMeshes, featureLevels);
-
-
     label iter = 0;
 
-    if (featureMeshes.size() && maxIter > 0)
+    if (meshRefiner_.features().size() && maxIter > 0)
     {
         for (; iter < maxIter; iter++)
         {
@@ -160,9 +96,6 @@ Foam::label Foam::autoRefineDriver::featureEdgeRefine
                     refineParams.keepPoints()[0],    // For now only use one.
                     refineParams.curvature(),
 
-                    featureMeshes,
-                    featureLevels,
-
                     true,               // featureRefinement
                     false,              // internalRefinement
                     false,              // surfaceRefinement
@@ -266,8 +199,6 @@ Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
         // Only look at surface intersections (minLevel and surface curvature),
         // do not do internal refinement (refinementShells)
 
-        const PtrList<featureEdgeMesh> dummyFeatures;
-
         labelList candidateCells
         (
             meshRefiner_.refineCandidates
@@ -275,9 +206,6 @@ Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
                 refineParams.keepPoints()[0],
                 refineParams.curvature(),
 
-                dummyFeatures,      // dummy featureMeshes;
-                labelList(0),       // dummy featureLevels;
-
                 false,              // featureRefinement
                 false,              // internalRefinement
                 true,               // surfaceRefinement
@@ -432,8 +360,6 @@ Foam::label Foam::autoRefineDriver::shellRefine
             << "----------------------------" << nl
             << endl;
 
-        const PtrList<featureEdgeMesh> dummyFeatures;
-
         labelList candidateCells
         (
             meshRefiner_.refineCandidates
@@ -441,9 +367,6 @@ Foam::label Foam::autoRefineDriver::shellRefine
                 refineParams.keepPoints()[0],
                 refineParams.curvature(),
 
-                dummyFeatures,      // dummy featureMeshes;
-                labelList(0),       // dummy featureLevels;
-
                 false,              // featureRefinement
                 true,               // internalRefinement
                 false,              // surfaceRefinement
@@ -459,7 +382,7 @@ Foam::label Foam::autoRefineDriver::shellRefine
                 << " cells to cellSet candidateCellsFromShells." << endl;
 
             cellSet c(mesh, "candidateCellsFromShells", candidateCells);
-            c.instance() = mesh.time().timeName();
+            c.instance() = meshRefiner_.timeName();
             c.write();
         }
 
@@ -677,6 +600,7 @@ void Foam::autoRefineDriver::splitAndMergeBaffles
         handleSnapProblems,                 // remove perp edge connected cells
         perpAngle,                          // perp angle
         false,                              // merge free standing baffles?
+        //true,                               // merge free standing baffles?
         motionDict,
         const_cast<Time&>(mesh.time()),
         globalToPatch_,
@@ -735,7 +659,8 @@ void Foam::autoRefineDriver::splitAndMergeBaffles
 
 void Foam::autoRefineDriver::mergePatchFaces
 (
-    const refinementParameters& refineParams
+    const refinementParameters& refineParams,
+    const dictionary& motionDict
 )
 {
     const fvMesh& mesh = meshRefiner_.mesh();
@@ -750,11 +675,12 @@ void Foam::autoRefineDriver::mergePatchFaces
         const_cast<Time&>(mesh.time())++;
     }
 
-    meshRefiner_.mergePatchFaces
+    meshRefiner_.mergePatchFacesUndo
     (
         Foam::cos(degToRad(45.0)),
         Foam::cos(degToRad(45.0)),
-        meshRefiner_.meshedPatches()
+        meshRefiner_.meshedPatches(),
+        motionDict
     );
 
     if (debug)
@@ -762,7 +688,7 @@ void Foam::autoRefineDriver::mergePatchFaces
         meshRefiner_.checkData();
     }
 
-    meshRefiner_.mergeEdges(Foam::cos(degToRad(45.0)));
+    meshRefiner_.mergeEdgesUndo(Foam::cos(degToRad(45.0)), motionDict);
 
     if (debug)
     {
@@ -789,13 +715,10 @@ void Foam::autoRefineDriver::doRefine
     // Check that all the keep points are inside the mesh.
     refineParams.findCells(mesh);
 
-    PtrList<dictionary> featDicts(refineDict.lookup("features"));
-
     // Refine around feature edges
     featureEdgeRefine
     (
         refineParams,
-        featDicts,
         100,    // maxIter
         0       // min cells to refine
     );
@@ -833,7 +756,7 @@ void Foam::autoRefineDriver::doRefine
     // Do something about cells with refined faces on the boundary
     if (prepareForSnapping)
     {
-        mergePatchFaces(refineParams);
+        mergePatchFaces(refineParams, motionDict);
     }
 
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -42,7 +42,6 @@ namespace Foam
 {
 
 // Forward declaration of classes
-class featureEdgeMesh;
 class refinementParameters;
 class meshRefinement;
 class decompositionMethod;
@@ -71,19 +70,10 @@ class autoRefineDriver
 
     // Private Member Functions
 
-        //- Read explicit feature edges
-        label readFeatureEdges
-        (
-            const PtrList<dictionary>& featDicts,
-            PtrList<featureEdgeMesh>& featureMeshes,
-            labelList& featureLevel
-        ) const;
-
         //- Refine all cells pierced by explicit feature edges
         label featureEdgeRefine
         (
             const refinementParameters& refineParams,
-            const PtrList<dictionary>& featDicts,
             const label maxIter,
             const label minRefine
         );
@@ -130,7 +120,8 @@ class autoRefineDriver
         //- Merge refined boundary faces (from exposing coarser cell)
         void mergePatchFaces
         (
-            const refinementParameters& refineParams
+            const refinementParameters& refineParams,
+            const dictionary& motionDict
         );
 
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,7 +27,6 @@ Description
 \*----------------------------------------------------------------------------*/
 
 #include "autoSnapDriver.H"
-#include "Time.H"
 #include "motionSmoother.H"
 #include "polyTopoChange.H"
 #include "OFstream.H"
@@ -36,13 +35,12 @@ Description
 #include "Time.H"
 #include "OFstream.H"
 #include "mapPolyMesh.H"
-#include "motionSmoother.H"
 #include "pointEdgePoint.H"
 #include "PointEdgeWave.H"
 #include "mergePoints.H"
 #include "snapParameters.H"
 #include "refinementSurfaces.H"
-
+#include "unitConversion.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -509,9 +507,7 @@ void Foam::autoSnapDriver::dumpMove
 )
 {
     // Dump direction of growth into file
-    Pout<< nl << "Dumping move direction to " << fName << nl
-        << "View this Lightwave-OBJ file with e.g. javaview" << nl
-        << endl;
+    Pout<< nl << "Dumping move direction to " << fName << endl;
 
     OFstream nearestStream(fName);
 
@@ -652,8 +648,7 @@ Foam::scalarField Foam::autoSnapDriver::calcSnapDistance
         -GREAT              // null value
     );
 
-    tmp<scalarField> tfld = snapParams.snapTol()*maxEdgeLen;
-    return tfld();
+    return snapParams.snapTol()*maxEdgeLen;
 }
 
 
@@ -688,6 +683,7 @@ void Foam::autoSnapDriver::preSmoothPatch
 
         // The current mesh is the starting mesh to smooth from.
         meshMover.setDisplacement(patchDisp);
+
         meshMover.correct();
 
         scalar oldErrorReduction = -1;
@@ -828,7 +824,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
             (
                 unzonedSurfaces,
                 localPoints,
-                sqr(4*snapDist),        // sqr of attract distance
+                sqr(snapDist),        // sqr of attract distance
                 hitSurface,
                 hitInfo
             );
@@ -880,7 +876,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
             (
                 labelList(1, zoneSurfI),
                 pointField(localPoints, zonePointIndices),
-                sqr(4*scalarField(minSnapDist, zonePointIndices)),
+                sqr(scalarField(minSnapDist, zonePointIndices)),
                 hitSurface,
                 hitInfo
             );
@@ -956,8 +952,8 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
         mesh,
         pp.meshPoints(),
         patchDisp,
-        minMagEqOp(),                   // combine op
-        vector(GREAT, GREAT, GREAT)     // null value
+        minMagSqrEqOp<point>(),         // combine op
+        vector(GREAT, GREAT, GREAT)     // null value (note: cannot use VGREAT)
     );
 
 
@@ -990,13 +986,12 @@ void Foam::autoSnapDriver::smoothDisplacement
 ) const
 {
     const fvMesh& mesh = meshRefiner_.mesh();
-    const pointMesh& pMesh = meshMover.pMesh();
     const indirectPrimitivePatch& pp = meshMover.patch();
 
     Info<< "Smoothing displacement ..." << endl;
 
     // Set edge diffusivity as inverse of distance to patch
-    scalarField edgeGamma(1.0/(edgePatchDist(pMesh, pp) + SMALL));
+    scalarField edgeGamma(1.0/(edgePatchDist(meshMover.pMesh(), pp) + SMALL));
     //scalarField edgeGamma(mesh.nEdges(), 1.0);
     //scalarField edgeGamma(wallGamma(mesh, pp, 10, 1));
 
@@ -1010,7 +1005,6 @@ void Foam::autoSnapDriver::smoothDisplacement
             Info<< "Iteration " << iter << endl;
         }
         pointVectorField oldDisp(disp);
-
         meshMover.smooth(oldDisp, edgeGamma, disp);
     }
     Info<< "Displacement smoothed in = "
@@ -1045,7 +1039,7 @@ void Foam::autoSnapDriver::smoothDisplacement
 }
 
 
-void Foam::autoSnapDriver::scaleMesh
+bool Foam::autoSnapDriver::scaleMesh
 (
     const snapParameters& snapParams,
     const label nInitErrors,
@@ -1061,6 +1055,8 @@ void Foam::autoSnapDriver::scaleMesh
 
     scalar oldErrorReduction = -1;
 
+    bool meshOk = false;
+
     Info<< "Moving mesh ..." << endl;
     for (label iter = 0; iter < 2*snapParams.nSnap(); iter++)
     {
@@ -1072,10 +1068,11 @@ void Foam::autoSnapDriver::scaleMesh
             oldErrorReduction = meshMover.setErrorReduction(0.0);
         }
 
-        if (meshMover.scaleMesh(checkFaces, baffles, true, nInitErrors))
+        meshOk = meshMover.scaleMesh(checkFaces, baffles, true, nInitErrors);
+
+        if (meshOk)
         {
             Info<< "Successfully moved mesh" << endl;
-
             break;
         }
         if (debug)
@@ -1098,6 +1095,8 @@ void Foam::autoSnapDriver::scaleMesh
     }
     Info<< "Moved mesh in = "
         << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
+
+    return meshOk;
 }
 
 
@@ -1191,7 +1190,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::repatchToSurface
         (
             unzonedSurfaces,
             localFaceCentres,
-            sqr(4*faceSnapDist),    // sqr of attract distance
+            sqr(faceSnapDist),    // sqr of attract distance
             hitSurface,
             hitRegion
         );
@@ -1260,6 +1259,7 @@ void Foam::autoSnapDriver::doSnap
 (
     const dictionary& snapDict,
     const dictionary& motionDict,
+    const scalar featureCos,
     const snapParameters& snapParams
 )
 {
@@ -1278,6 +1278,21 @@ void Foam::autoSnapDriver::doSnap
     List<labelPair> baffles;
     meshRefiner_.createZoneBaffles(globalToPatch_, baffles);
 
+
+    bool doFeatures = false;
+    label nFeatIter = 1;
+    if (snapParams.nFeatureSnap() > 0)
+    {
+        doFeatures = true;
+        nFeatIter = snapParams.nFeatureSnap();
+
+        Info<< "Snapping to features in " << nFeatIter
+            << " iterations ..." << endl;
+    }
+
+
+    bool meshOk = false;
+
     {
         autoPtr<indirectPrimitivePatch> ppPtr
         (
@@ -1287,10 +1302,9 @@ void Foam::autoSnapDriver::doSnap
                 adaptPatchIDs
             )
         );
-        indirectPrimitivePatch& pp = ppPtr();
 
         // Distance to attract to nearest feature on surface
-        const scalarField snapDist(calcSnapDistance(snapParams, pp));
+        const scalarField snapDist(calcSnapDistance(snapParams, ppPtr()));
 
 
         // Construct iterative mesh mover.
@@ -1302,7 +1316,7 @@ void Foam::autoSnapDriver::doSnap
         motionSmoother meshMover
         (
             mesh,
-            pp,
+            ppPtr(),
             adaptPatchIDs,
             meshRefinement::makeDisplacementField(pMesh, adaptPatchIDs),
             motionDict
@@ -1330,19 +1344,97 @@ void Foam::autoSnapDriver::doSnap
         // Pre-smooth patch vertices (so before determining nearest)
         preSmoothPatch(snapParams, nInitErrors, baffles, meshMover);
 
-        // Calculate displacement at every patch point. Insert into
-        // meshMover.
-        calcNearestSurface(snapDist, meshMover);
+        if (debug)
+        {
+            Pout<< "Writing patch smoothed mesh to time "
+                << meshRefiner_.timeName() << '.' << endl;
+            meshRefiner_.write
+            (
+                debug,
+                mesh.time().path()/meshRefiner_.timeName()
+            );
+            Pout<< "Dumped mesh in = "
+                << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
+        }
 
-        //// Get smoothly varying internal displacement field.
-        //- 2009-12-16 : was not found to be beneficial. Keeping internal
-        // fields fixed slightly increases skewness (on boundary)
-        // but lowers non-orthogonality quite a bit (e.g. 65->59 degrees).
-        // Maybe if better smoother?
-        //smoothDisplacement(snapParams, meshMover);
 
-        // Apply internal displacement to mesh.
-        scaleMesh(snapParams, nInitErrors, baffles, meshMover);
+        for (label iter = 0; iter < nFeatIter; iter++)
+        {
+            Info<< nl
+                << "Morph iteration " << iter << nl
+                << "-----------------" << endl;
+
+            // Calculate displacement at every patch point. Insert into
+            // meshMover.
+            vectorField disp = calcNearestSurface(snapDist, meshMover);
+
+            // Override displacement with feature edge attempt
+            if (doFeatures)
+            {
+                disp = calcNearestSurfaceFeature
+                (
+                    iter,
+                    featureCos,
+                    scalar(iter+1)/nFeatIter,
+                    snapDist,
+                    disp,
+                    meshMover
+                );
+            }
+
+            if (debug&meshRefinement::OBJINTERSECTIONS)
+            {
+                dumpMove
+                (
+                    mesh.time().path()
+                  / "patchDisplacement_" + name(iter) + ".obj",
+                    ppPtr().localPoints(),
+                    ppPtr().localPoints() + disp
+                );
+            }
+
+            // Get smoothly varying internal displacement field.
+            smoothDisplacement(snapParams, meshMover);
+
+            // Apply internal displacement to mesh.
+            meshOk = scaleMesh
+            (
+                snapParams,
+                nInitErrors,
+                baffles,
+                meshMover
+            );
+
+            if (!meshOk)
+            {
+                Info<< "Did not succesfully snap mesh. Giving up."
+                    << nl << endl;
+
+                // Use current mesh as base mesh
+                meshMover.correct();
+
+                break;
+            }
+
+            if (debug)
+            {
+                const_cast<Time&>(mesh.time())++;
+                Pout<< "Writing scaled mesh to time "
+                    << meshRefiner_.timeName() << endl;
+                meshRefiner_.write
+                (
+                    debug,
+                    mesh.time().path()/meshRefiner_.timeName()
+                );
+                Pout<< "Writing displacement field ..." << endl;
+                meshMover.displacement().write();
+                tmp<pointScalarField> magDisp(mag(meshMover.displacement()));
+                magDisp().write();
+            }
+
+            // Use current mesh as base mesh
+            meshMover.correct();
+        }
     }
 
     // Merge any introduced baffles.
@@ -1350,6 +1442,34 @@ void Foam::autoSnapDriver::doSnap
 
     // Repatch faces according to nearest.
     repatchToSurface(snapParams, adaptPatchIDs);
+
+    // Repatching might have caused faces to be on same patch and hence
+    // mergeable so try again to merge coplanar faces
+    label nChanged = meshRefiner_.mergePatchFacesUndo
+    (
+        featureCos,  // minCos
+        featureCos,  // concaveCos
+        meshRefiner_.meshedPatches(),
+        motionDict
+    );
+
+    nChanged += meshRefiner_.mergeEdgesUndo
+    (
+        featureCos,
+        motionDict
+    );
+
+    if (nChanged > 0 && debug)
+    {
+        const_cast<Time&>(mesh.time())++;
+        Pout<< "Writing patchFace merged mesh to time "
+            << meshRefiner_.timeName() << endl;
+        meshRefiner_.write
+        (
+            debug,
+            mesh.time().path()/meshRefiner_.timeName()
+        );
+    }
 }
 
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,6 +29,7 @@ Description
 
 SourceFiles
     autoSnapDriver.C
+    autoSnapDriverFeature.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -45,6 +46,7 @@ namespace Foam
 // Forward declaration of classes
 class motionSmoother;
 class snapParameters;
+class pointConstraint;
 
 /*---------------------------------------------------------------------------*\
                            Class autoSnapDriver Declaration
@@ -52,23 +54,6 @@ class snapParameters;
 
 class autoSnapDriver
 {
-    // Private classes
-
-        //- Combine operator class for equalizing displacements.
-        class minMagEqOp
-        {
-        public:
-
-            void operator()(vector& x, const vector& y) const
-            {
-                if (magSqr(y) < magSqr(x))
-                {
-                    x = y;
-                }
-            }
-        };
-
-
     // Private data
 
         //- Mesh+surface
@@ -124,6 +109,174 @@ class autoSnapDriver
                 const vectorField&
             );
 
+            // Feature line snapping
+
+                void smoothAndConstrain
+                (
+                    const indirectPrimitivePatch& pp,
+                    const List<pointConstraint>& constraints,
+                    vectorField& disp
+                ) const;
+                void calcNearest
+                (
+                    const pointField& points,
+                    vectorField& disp,
+                    vectorField& surfaceNormal
+                ) const;
+                void calcNearestFace
+                (
+                    const label iter,
+                    const indirectPrimitivePatch& pp,
+                    vectorField& faceDisp,
+                    vectorField& faceSurfaceNormal,
+                    vectorField& faceRotation
+                ) const;
+                void interpolateFaceToPoint
+                (
+                    const label iter,
+                    const indirectPrimitivePatch& pp,
+                    const vectorField& faceSurfaceNormal,
+
+                    const vectorField& faceDisp,
+                    const vectorField& faceRotation,
+
+                    vectorField& patchDisp,
+                    vectorField& patchRotationDisp
+                ) const;
+                void correctAttraction
+                (
+                    const DynamicList<point>& surfacePoints,
+                    const DynamicList<label>& surfaceCount,
+                    const point& edgePt,
+                    const vector& edgeNormal,   // normalised normal
+                    const point& pt,
+                    vector& edgeOffset  // offset from pt to point on edge
+                ) const;
+                void binFeatureFace
+                (
+                    const label iter,
+                    const scalar featureCos,
+
+                    const indirectPrimitivePatch& pp,
+                    const scalar snapDist,
+
+                    const point& fc,
+                    const vector& faceSurfaceNormal,
+                    const vector& faceDisp,
+
+                    DynamicList<point>& surfacePoints,
+                    DynamicList<vector>& surfaceNormals,
+                    DynamicList<label>& surfaceCount
+                ) const;
+                void binFeatureFaces
+                (
+                    const label iter,
+                    const scalar featureCos,
+
+                    const indirectPrimitivePatch& pp,
+                    const scalarField& snapDist,
+                    const label pointI,
+
+                    const List<List<point> >& pointFaceNormals,
+                    const List<List<point> >& pointFaceDisp,
+                    const List<List<point> >& pointFaceCentres,
+
+                    DynamicList<point>& surfacePoints,
+                    DynamicList<vector>& surfaceNormals,
+                    DynamicList<label>& surfaceCount
+                ) const;
+
+                void featureAttractionUsingReconstruction
+                (
+                    const label iter,
+                    const scalar featureCos,
+                    const indirectPrimitivePatch& pp,
+                    const scalarField& snapDist,
+
+                    const List<List<point> >& pointFaceNormals,
+                    const List<List<point> >& pointFaceDisp,
+                    const List<List<point> >& pointFaceCentres,
+
+                    vectorField& patchAttraction,
+                    List<pointConstraint>& patchConstraints
+                ) const;
+
+                void determineAllFeatures
+                (
+                    const label iter,
+                    const scalar featureCos,
+
+                    const indirectPrimitivePatch&,
+                    const scalarField&,
+
+                    const List<List<point> >& pointFaceNormals,
+                    const List<List<point> >& pointFaceDisp,
+                    const List<List<point> >& pointFaceCentres,
+
+                    List<labelList>& pointAttractor,
+                    List<List<pointConstraint> >& pointConstraints,
+                    // Feature-edge to pp point
+                    List<List<DynamicList<point> > >& edgeAttractors,
+                    List<List<DynamicList<pointConstraint> > >& edgeConstraints,
+                    vectorField& patchAttraction,
+                    List<pointConstraint>& patchConstraints
+                ) const;
+                void determineFeatures
+                (
+                    const label iter,
+                    const scalar featureCos,
+
+                    const indirectPrimitivePatch&,
+                    const scalarField&,
+
+                    const List<List<point> >& pointFaceNormals,
+                    const List<List<point> >& pointFaceDisp,
+                    const List<List<point> >& pointFaceCentres,
+
+                    List<labelList>& pointAttractor,
+                    List<List<pointConstraint> >& pointConstraints,
+                    // Feature-edge to pp point
+                    List<List<DynamicList<point> > >& edgeAttractors,
+                    List<List<DynamicList<pointConstraint> > >& edgeConstraints,
+                    vectorField& patchAttraction,
+                    List<pointConstraint>& patchConstraints
+                ) const;
+
+                void featureAttractionUsingFeatureEdges
+                (
+                    const label iter,
+                    const scalar featureCos,
+                    const indirectPrimitivePatch& pp,
+                    const scalarField& snapDist,
+
+                    const List<List<point> >& pointFaceNormals,
+                    const List<List<point> >& pointFaceDisp,
+                    const List<List<point> >& pointFaceCentres,
+
+                    vectorField& patchAttraction,
+                    List<pointConstraint>& patchConstraints
+                ) const;
+                void preventFaceSqueeze
+                (
+                    const label iter,
+                    const scalar featureCos,
+                    const indirectPrimitivePatch& pp,
+                    const scalarField& snapDist,
+
+                    vectorField& patchAttraction,
+                    List<pointConstraint>& patchConstraints
+                ) const;
+
+                vectorField calcNearestSurfaceFeature
+                (
+                    const label iter,
+                    const scalar featureCos,
+                    const scalar featureAttract,
+                    const scalarField& snapDist,
+                    const vectorField& nearestDisp,
+                    motionSmoother& meshMover
+                ) const;
+
 
         //- Disallow default bitwise copy construct
         autoSnapDriver(const autoSnapDriver&);
@@ -188,6 +341,7 @@ public:
                 motionSmoother& meshMover
             ) const;
 
+
             //- Smooth the displacement field to the internal.
             void smoothDisplacement
             (
@@ -196,8 +350,9 @@ public:
             ) const;
 
             //- Do the hard work: move the mesh according to displacement,
-            //  locally relax the displacement.
-            void scaleMesh
+            //  locally relax the displacement. Return true if ended up with
+            //  correct mesh, false if not.
+            bool scaleMesh
             (
                 const snapParameters& snapParams,
                 const label nInitErrors,
@@ -212,11 +367,11 @@ public:
                 const labelList& adaptPatchIDs
             );
 
-
             void doSnap
             (
                 const dictionary& snapDict,
                 const dictionary& motionDict,
+                const scalar featureCos,
                 const snapParameters& snapParams
             );
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,23 +27,14 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from dictionary
-Foam::snapParameters::snapParameters(const dictionary& dict, const label dummy)
-:
-    nSmoothPatch_(readLabel(dict.lookup("nSmoothPatch"))),
-    snapTol_(readScalar(dict.lookup("snapTol"))),
-    nSmoothDispl_(readLabel(dict.lookup("nSmoothDispl"))),
-    nSnap_(readLabel(dict.lookup("nSnap")))
-{}
-
-
 // Construct from dictionary
 Foam::snapParameters::snapParameters(const dictionary& dict)
 :
     nSmoothPatch_(readLabel(dict.lookup("nSmoothPatch"))),
     snapTol_(readScalar(dict.lookup("tolerance"))),
     nSmoothDispl_(readLabel(dict.lookup("nSolveIter"))),
-    nSnap_(readLabel(dict.lookup("nRelaxIter")))
+    nSnap_(readLabel(dict.lookup("nRelaxIter"))),
+    nFeatureSnap_(dict.lookupOrDefault("nFeatureSnapIter", -1))
 {}
 
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -61,6 +61,8 @@ class snapParameters
 
         const label nSnap_;
 
+        const label nFeatureSnap_;
+
 
     // Private Member Functions
 
@@ -75,10 +77,7 @@ public:
 
     // Constructors
 
-        //- Construct from dictionary - old syntax
-        snapParameters(const dictionary& dict, const label dummy);
-
-        //- Construct from dictionary - new syntax
+        //- Construct from dictionary
         snapParameters(const dictionary& dict);
 
 
@@ -114,6 +113,12 @@ public:
             {
                 return nSnap_;
             }
+
+            label nFeatureSnap() const
+            {
+                return nFeatureSnap_;
+            }
+
 };
 
 

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

@@ -31,6 +31,7 @@ License
 #include "Time.H"
 #include "refinementHistory.H"
 #include "refinementSurfaces.H"
+#include "refinementFeatures.H"
 #include "decompositionMethod.H"
 #include "regionSplit.H"
 #include "fvMeshDistribute.H"
@@ -497,6 +498,10 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::doRemoveCells
         mesh_.clearOut();
     }
 
+    // Reset the instance for if in overwrite mode
+    mesh_.setInstance(timeName());
+    setInstance(mesh_.facesInstance());
+
     // Update local mesh data
     cellRemover.updateMesh(map);
 
@@ -857,6 +862,7 @@ Foam::meshRefinement::meshRefinement
     const scalar mergeDistance,
     const bool overwrite,
     const refinementSurfaces& surfaces,
+    const refinementFeatures& features,
     const shellSurfaces& shells
 )
 :
@@ -865,6 +871,7 @@ Foam::meshRefinement::meshRefinement
     overwrite_(overwrite),
     oldInstance_(mesh.pointsInstance()),
     surfaces_(surfaces),
+    features_(features),
     shells_(shells),
     meshCutter_
     (
@@ -1271,6 +1278,10 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
 
         // Update numbering of meshRefiner
         distribute(map);
+
+        // Set correct instance (for if overwrite)
+        mesh_.setInstance(timeName());
+        setInstance(mesh_.facesInstance());
     }
     return map;
 }
@@ -1449,6 +1460,8 @@ Foam::tmp<Foam::pointVectorField> Foam::meshRefinement::makeDisplacementField
         }
     }
 
+    // Note: time().timeName() instead of meshRefinement::timeName() since
+    // postprocessable field.
     tmp<pointVectorField> tfld
     (
         new pointVectorField
@@ -1456,7 +1469,7 @@ Foam::tmp<Foam::pointVectorField> Foam::meshRefinement::makeDisplacementField
             IOobject
             (
                 "pointDisplacement",
-                mesh.time().timeName(), //timeName(),
+                mesh.time().timeName(),
                 mesh,
                 IOobject::NO_READ,
                 IOobject::AUTO_WRITE
@@ -1930,10 +1943,6 @@ void Foam::meshRefinement::distribute(const mapDistributePolyMesh& map)
             pointMap.clear();
         }
     }
-
-    // If necessary reset the instance
-    mesh_.setInstance(timeName());
-    setInstance(mesh_.facesInstance());
 }
 
 
@@ -2059,10 +2068,6 @@ void Foam::meshRefinement::updateMesh
             data.transfer(newFaceData);
         }
     }
-
-    // If necessary reset the instance
-    mesh_.setInstance(timeName());
-    setInstance(mesh_.facesInstance());
 }
 
 
@@ -2118,12 +2123,32 @@ void Foam::meshRefinement::printMeshInfo(const bool debug, const string& msg)
             << "  points(local):" << mesh_.nPoints()
             << endl;
     }
-    else
+
     {
+        PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
+        label nMasterFaces = 0;
+        forAll(isMasterFace, i)
+        {
+            if (isMasterFace[i])
+            {
+                nMasterFaces++;
+            }
+        }
+
+        PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh_));
+        label nMasterPoints = 0;
+        forAll(isMasterPoint, i)
+        {
+            if (isMasterPoint[i])
+            {
+                nMasterPoints++;
+            }
+        }
+
         Info<< msg.c_str()
             << " : cells:" << pData.nTotalCells()
-            << "  faces:" << pData.nTotalFaces()
-            << "  points:" << pData.nTotalPoints()
+            << "  faces:" << returnReduce(nMasterFaces, sumOp<label>())
+            << "  points:" << returnReduce(nMasterPoints, sumOp<label>())
             << endl;
     }
 
@@ -2168,12 +2193,14 @@ Foam::word Foam::meshRefinement::timeName() const
 
 void Foam::meshRefinement::dumpRefinementLevel() const
 {
+    // Note: use time().timeName(), not meshRefinement::timeName()
+    // so as to dump the fields to 0, not to constant.
     volScalarField volRefLevel
     (
         IOobject
         (
             "cellLevel",
-            mesh_.time().timeName(),// Dump to current time, not to mesh inst
+            mesh_.time().timeName(),
             mesh_,
             IOobject::NO_READ,
             IOobject::AUTO_WRITE,
@@ -2201,7 +2228,7 @@ void Foam::meshRefinement::dumpRefinementLevel() const
         IOobject
         (
             "pointLevel",
-            mesh_.time().timeName(),// Dump to current time, not to mesh inst
+            mesh_.time().timeName(),
             mesh_,
             IOobject::NO_READ,
             IOobject::NO_WRITE,

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

@@ -62,13 +62,14 @@ class fvMesh;
 class mapDistributePolyMesh;
 class decompositionMethod;
 class refinementSurfaces;
+class refinementFeatures;
 class shellSurfaces;
 class removeCells;
-class featureEdgeMesh;
 class fvMeshDistribute;
 class searchableSurface;
 class regionSplit;
 class globalIndex;
+class removePoints;
 
 /*---------------------------------------------------------------------------*\
                            Class meshRefinement Declaration
@@ -118,6 +119,9 @@ private:
         //- All surface-intersection interaction
         const refinementSurfaces& surfaces_;
 
+        //- All feature-edge interaction
+        const refinementFeatures& features_;
+
         //- All shell-refinement interaction
         const shellSurfaces& shells_;
 
@@ -260,8 +264,6 @@ private:
             label markFeatureRefinement
             (
                 const point& keepPoint,
-                const PtrList<featureEdgeMesh>& featureMeshes,
-                const labelList& featureLevels,
                 const label nAllowRefine,
 
                 labelList& refineCell,
@@ -415,12 +417,6 @@ private:
                 const labelList& globalToPatch
             ) const;
 
-            ////- Initial test of marking faces using geometric information.
-            //labelList markFacesOnProblemCellsGeometric
-            //(
-            //    const dictionary& motionDict
-            //) const;
-
 
         // Baffle merging
 
@@ -501,6 +497,7 @@ public:
             const scalar mergeDistance,
             const bool overwrite,
             const refinementSurfaces&,
+            const refinementFeatures&,
             const shellSurfaces&
         );
 
@@ -542,6 +539,12 @@ public:
                 return surfaces_;
             }
 
+            //- reference to feature edge mesh
+            const refinementFeatures& features() const
+            {
+                return features_;
+            }
+
             //- reference to refinement shells (regions)
             const shellSurfaces& shells() const
             {
@@ -648,9 +651,6 @@ public:
                 const point& keepPoint,
                 const scalar curvature,
 
-                const PtrList<featureEdgeMesh>& featureMeshes,
-                const labelList& featureLevels,
-
                 const bool featureRefinement,
                 const bool internalRefinement,
                 const bool surfaceRefinement,
@@ -803,19 +803,45 @@ public:
                     const Map<label>& cellsToRestore
                 );
 
+            // Merging coplanar faces and edges
 
-            //- Merge faces on the same patch (usually from exposing refinement)
-            //  Returns global number of faces merged.
-            label mergePatchFaces
-            (
-                const scalar minCos,
-                const scalar concaveCos,
-                const labelList& patchIDs
-            );
+                label mergePatchFacesUndo
+                (
+                    const scalar minCos,
+                    const scalar concaveCos,
+                    const labelList& patchIDs,
+                    const dictionary& motionDict
+                );
 
-            //- Remove points not used by any face or points used
-            // by only two faces where the edges are in line
-            autoPtr<mapPolyMesh> mergeEdges(const scalar minCos);
+                autoPtr<mapPolyMesh> doRemovePoints
+                (
+                    removePoints& pointRemover,
+                    const boolList& pointCanBeDeleted
+                );
+
+                autoPtr<mapPolyMesh> doRestorePoints
+                (
+                    removePoints& pointRemover,
+                    const labelList& facesToRestore
+                );
+
+                labelList collectFaces
+                (
+                    const labelList& candidateFaces,
+                    const labelHashSet& set
+                ) const;
+
+                // Pick up faces of cells of faces in set.
+                labelList growFaceCellFace
+                (
+                    const labelHashSet& set
+                ) const;
+
+                label mergeEdgesUndo
+                (
+                    const scalar minCos,
+                    const dictionary& motionDict
+                );
 
 
         // Debug/IO

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

@@ -524,6 +524,10 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles
         mesh_.clearOut();
     }
 
+
+    // Reset the instance for if in overwrite mode
+    mesh_.setInstance(timeName());
+
     //- Redo the intersections on the newly create baffle faces. Note that
     //  this changes also the cell centre positions.
     faceSet baffledFacesSet(mesh_, "baffledFacesSet", 2*nBaffles);
@@ -644,7 +648,7 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::getDuplicateFaces
         }
         Pout<< "Writing duplicate faces (baffles) to faceSet "
             << duplicateFaceSet.name() << nl << endl;
-        duplicateFaceSet.instance() = mesh_.time().timeName();
+        duplicateFaceSet.instance() = timeName();
         duplicateFaceSet.write();
     }
 
@@ -724,9 +728,9 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createZoneBaffles
             if (debug)
             {
                 const_cast<Time&>(mesh_.time())++;
-                Pout<< "Writing baffled mesh to time "
-                    << mesh_.time().timeName() << endl;
-                mesh_.write();
+                Pout<< "Writing zone-baffled mesh to time " << timeName()
+                    << endl;
+                write(debug, mesh_.time().path()/"baffles");
             }
         }
         Info<< "Created " << nZoneFaces << " baffles in = "
@@ -746,14 +750,6 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::filterDuplicateFaces
     const List<labelPair>& couples
 ) const
 {
-    // Construct addressing engine for all duplicate faces (only one
-    // for each pair)
-
-    // Duplicate faces in mesh labels (first face of each pair only)
-    // (reused later on to mark off filtered couples. see below)
-    labelList duplicateFaces(couples.size());
-
-
     // All duplicate faces on edge of the patch are to be merged.
     // So we count for all edges of duplicate faces how many duplicate
     // faces use them.
@@ -824,7 +820,8 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::filterDuplicateFaces
 
 
     // Baffles which are not next to other boundaries and baffles will have
-    // value 2*1000000+2*1
+    // nBafflesPerEdge value 2*1000000+2*1 (from 2 boundary faces which are
+    // both baffle faces)
 
     List<labelPair> filteredCouples(couples.size());
     label filterI = 0;
@@ -967,6 +964,9 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::mergeBaffles
         mesh_.clearOut();
     }
 
+    // Reset the instance for if in overwrite mode
+    mesh_.setInstance(timeName());
+
     // Update intersections. Recalculate intersections on merged faces since
     // this seems to give problems? Note: should not be nessecary since
     // baffles preserve intersections from when they were created.
@@ -1789,7 +1789,7 @@ void Foam::meshRefinement::baffleAndSplitMesh
                     problemTopo.insert(faceI);
                 }
             }
-            problemTopo.instance() = mesh_.time().timeName();
+            problemTopo.instance() = timeName();
             Pout<< "Dumping " << problemTopo.size()
                 << " problem faces to " << problemTopo.objectPath() << endl;
             problemTopo.write();
@@ -2246,6 +2246,9 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::dupNonManifoldPoints()
         mesh_.clearOut();
     }
 
+    // Reset the instance for if in overwrite mode
+    mesh_.setInstance(timeName());
+
     // Update intersections. Is mapping only (no faces created, positions stay
     // same) so no need to recalculate intersections.
     updateMesh(map, labelList(0));
@@ -2782,6 +2785,9 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
         mesh_.clearOut();
     }
 
+    // Reset the instance for if in overwrite mode
+    mesh_.setInstance(timeName());
+
     // Print some stats (note: zones are synchronised)
     if (mesh_.cellZones().size() > 0)
     {

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,42 +27,258 @@ License
 #include "combineFaces.H"
 #include "polyTopoChange.H"
 #include "removePoints.H"
+#include "faceSet.H"
+#include "Time.H"
+#include "motionSmoother.H"
+#include "syncTools.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-// Merge faces that are in-line.
-Foam::label Foam::meshRefinement::mergePatchFaces
+//// 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;
+//}
+//
+//
+//// Remove points not used by any face or points used by only two faces where
+//// the edges are in line
+//Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::mergeEdges
+//(
+//    const scalar minCos
+//)
+//{
+//    // Point removal analysis engine
+//    removePoints pointRemover(mesh_);
+//
+//    // Count usage of points
+//    boolList pointCanBeDeleted;
+//    label nRemove = pointRemover.countPointUsage(minCos, pointCanBeDeleted);
+//
+//    Info<< "Removing " << nRemove
+//        << " straight edge points." << endl;
+//
+//    autoPtr<mapPolyMesh> map;
+//
+//    if (nRemove > 0)
+//    {
+//        // Save my local faces that will change. These changed faces might
+//        // cause a shift in the cell centre which needs to be retested.
+//        // Have to do this before changing mesh since point will be removed.
+//        labelHashSet retestOldFaces(nRemove / Pstream::nProcs());
+//
+//        {
+//            const faceList& faces = mesh_.faces();
+//
+//            forAll(faces, faceI)
+//            {
+//                const face& f = faces[faceI];
+//
+//                forAll(f, fp)
+//                {
+//                    if (pointCanBeDeleted[f[fp]])
+//                    {
+//                        retestOldFaces.insert(faceI);
+//                        break;
+//                    }
+//                }
+//            }
+//        }
+//
+//        // Topology changes container
+//        polyTopoChange meshMod(mesh_);
+//
+//        pointRemover.setRefinement(pointCanBeDeleted, meshMod);
+//
+//        // Change the mesh (no inflation)
+//        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());
+//
+//        pointRemover.updateMesh(map);
+//
+//        // Get the kept faces that need to be recalculated.
+//        labelHashSet retestFaces(6*retestOldFaces.size());
+//
+//        const cellList& cells = mesh_.cells();
+//
+//        forAllConstIter(labelHashSet, retestOldFaces, iter)
+//        {
+//            label faceI = map().reverseFaceMap()[iter.key()];
+//
+//            const cell& ownFaces = cells[mesh_.faceOwner()[faceI]];
+//
+//            forAll(ownFaces, i)
+//            {
+//                retestFaces.insert(ownFaces[i]);
+//            }
+//
+//            if (mesh_.isInternalFace(faceI))
+//            {
+//                const cell& neiFaces = cells[mesh_.faceNeighbour()[faceI]];
+//
+//                forAll(neiFaces, i)
+//                {
+//                    retestFaces.insert(neiFaces[i]);
+//                }
+//            }
+//        }
+//        updateMesh(map, retestFaces.toc());
+//    }
+//
+//    return map;
+//}
+
+
+Foam::label Foam::meshRefinement::mergePatchFacesUndo
 (
     const scalar minCos,
     const scalar concaveCos,
-    const labelList& patchIDs
+    const labelList& patchIDs,
+    const dictionary& motionDict
 )
 {
     // Patch face merging engine
-    combineFaces faceCombiner(mesh_);
-
-    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+    combineFaces faceCombiner(mesh_, true);
 
     // Pick up all candidate cells on boundary
     labelHashSet boundaryCells(mesh_.nFaces()-mesh_.nInternalFaces());
 
-    forAll(patchIDs, i)
     {
-        label patchI = patchIDs[i];
+        const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
-        const polyPatch& patch = patches[patchI];
-
-        if (!patch.coupled())
+        forAll(patchIDs, i)
         {
-            forAll(patch, i)
+            label patchI = patchIDs[i];
+
+            const polyPatch& patch = patches[patchI];
+
+            if (!patch.coupled())
             {
-                boundaryCells.insert(mesh_.faceOwner()[patch.start()+i]);
+                forAll(patch, i)
+                {
+                    boundaryCells.insert(mesh_.faceOwner()[patch.start()+i]);
+                }
             }
         }
     }
 
     // Get all sets of faces that can be merged
-    labelListList mergeSets
+    labelListList allFaceSets
     (
         faceCombiner.getMergeSets
         (
@@ -72,19 +288,44 @@ Foam::label Foam::meshRefinement::mergePatchFaces
         )
     );
 
-    label nFaceSets = returnReduce(mergeSets.size(), sumOp<label>());
+    label nFaceSets = returnReduce(allFaceSets.size(), sumOp<label>());
 
-    Info<< "mergePatchFaces : Merging " << nFaceSets
-        << " sets of faces." << endl;
+    Info<< "Merging " << nFaceSets << " sets of faces." << nl << endl;
 
     if (nFaceSets > 0)
     {
+        if (debug)
+        {
+            faceSet allSets(mesh_, "allFaceSets", allFaceSets.size());
+            forAll(allFaceSets, setI)
+            {
+                forAll(allFaceSets[setI], i)
+                {
+                    allSets.insert(allFaceSets[setI][i]);
+                }
+            }
+            Pout<< "Writing all faces to be merged to set "
+                << allSets.objectPath() << endl;
+            allSets.instance() = timeName();
+            allSets.write();
+
+            const_cast<Time&>(mesh_.time())++;
+        }
+
+
         // 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);
+        // Merge all faces of a set into the first face of the set.
+        faceCombiner.setRefinement(allFaceSets, meshMod);
+
+        // Experimental: store data for all the points that have been deleted
+        storeData
+        (
+            faceCombiner.savedPointLabels(),    // points to store
+            labelList(0),                       // faces to store
+            labelList(0)                        // cells to store
+        );
 
         // Change the mesh (no inflation)
         autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
@@ -99,138 +340,538 @@ Foam::label Foam::meshRefinement::mergePatchFaces
         }
         else
         {
-            // Delete mesh volumes. No other way to do this?
+            // Delete mesh volumes.
             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());
+        updateMesh(map, labelList(0));
 
-        forAll(mergeSets, setI)
+        if (debug)
         {
-            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]);
-            }
+            Pout<< "Writing initial merged-faces mesh to time "
+                << timeName() << nl << endl;
+            write();
         }
-        updateMesh(map, retestFaces.toc());
-    }
 
+        for (label iteration = 0; iteration < 100; iteration++)
+        {
+            Info<< nl
+                << "Undo iteration " << iteration << nl
+                << "----------------" << endl;
+
+
+            // Check mesh for errors
+            // ~~~~~~~~~~~~~~~~~~~~~
+
+            faceSet errorFaces
+            (
+                mesh_,
+                "errorFaces",
+                mesh_.nFaces()-mesh_.nInternalFaces()
+            );
+            bool hasErrors = motionSmoother::checkMesh
+            (
+                false,  // report
+                mesh_,
+                motionDict,
+                errorFaces
+            );
+
+            //if (checkEdgeConnectivity)
+            //{
+            //    Info<< "Checking edge-face connectivity (duplicate faces"
+            //        << " or non-consecutive shared vertices)" << endl;
+            //
+            //    label nOldSize = errorFaces.size();
+            //
+            //    hasErrors =
+            //        mesh_.checkFaceFaces
+            //        (
+            //            false,
+            //            &errorFaces
+            //        )
+            //     || hasErrors;
+            //
+            //    Info<< "Detected additional "
+            //        <<  returnReduce
+            //            (
+            //                errorFaces.size() - nOldSize,
+            //                sumOp<label>()
+            //            )
+            //        << " faces with illegal face-face connectivity"
+            //        << endl;
+            //}
+
+            if (!hasErrors)
+            {
+                break;
+            }
+
+
+            if (debug)
+            {
+                errorFaces.instance() = timeName();
+                Pout<< "Writing all faces in error to faceSet "
+                    << errorFaces.objectPath() << nl << endl;
+                errorFaces.write();
+            }
+
+
+            // Check any master cells for using any of the error faces
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+            DynamicList<label> mastersToRestore(allFaceSets.size());
+
+            forAll(allFaceSets, setI)
+            {
+                label masterFaceI = faceCombiner.masterFace()[setI];
+
+                if (masterFaceI != -1)
+                {
+                    label masterCellII = mesh_.faceOwner()[masterFaceI];
+
+                    const cell& cFaces = mesh_.cells()[masterCellII];
+
+                    forAll(cFaces, i)
+                    {
+                        if (errorFaces.found(cFaces[i]))
+                        {
+                            mastersToRestore.append(masterFaceI);
+                            break;
+                        }
+                    }
+                }
+            }
+            mastersToRestore.shrink();
+
+            label nRestore = returnReduce
+            (
+                mastersToRestore.size(),
+                sumOp<label>()
+            );
+
+            Info<< "Masters that need to be restored:"
+                << nRestore << endl;
+
+            if (debug)
+            {
+                faceSet restoreSet(mesh_, "mastersToRestore", mastersToRestore);
+                restoreSet.instance() = timeName();
+                Pout<< "Writing all " << mastersToRestore.size()
+                    << " masterfaces to be restored to set "
+                    << restoreSet.objectPath() << endl;
+                restoreSet.write();
+            }
+
+
+            if (nRestore == 0)
+            {
+                break;
+            }
+
+
+            // Undo
+            // ~~~~
+
+            if (debug)
+            {
+                const_cast<Time&>(mesh_.time())++;
+            }
+
+            // Topology changes container
+            polyTopoChange meshMod(mesh_);
+
+            // Merge all faces of a set into the first face of the set.
+            // Experimental:mark all points/faces/cells that have been restored.
+            Map<label> restoredPoints(0);
+            Map<label> restoredFaces(0);
+            Map<label> restoredCells(0);
+
+            faceCombiner.setUnrefinement
+            (
+                mastersToRestore,
+                meshMod,
+                restoredPoints,
+                restoredFaces,
+                restoredCells
+            );
+
+            // Change the mesh (no inflation)
+            autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
+
+            // Update fields
+            mesh_.updateMesh(map);
+
+            // Move mesh (since morphing does not do this)
+            if (map().hasMotionPoints())
+            {
+                mesh_.movePoints(map().preMotionPoints());
+            }
+            else
+            {
+                // Delete mesh volumes.
+                mesh_.clearOut();
+            }
+
+
+            // Reset the instance for if in overwrite mode
+            mesh_.setInstance(timeName());
+
+            faceCombiner.updateMesh(map);
+
+            // Renumber restore maps
+            inplaceMapKey(map().reversePointMap(), restoredPoints);
+            inplaceMapKey(map().reverseFaceMap(), restoredFaces);
+            inplaceMapKey(map().reverseCellMap(), restoredCells);
+
+            // Experimental:restore all points/face/cells in maps
+            updateMesh
+            (
+                map,
+                labelList(0),       // changedFaces
+                restoredPoints,
+                restoredFaces,
+                restoredCells
+            );
+
+            if (debug)
+            {
+                Pout<< "Writing merged-faces mesh to time "
+                    << timeName() << nl << endl;
+                write();
+            }
+
+            Info<< endl;
+        }
+    }
+    else
+    {
+        Info<< "No faces merged ..." << endl;
+    }
 
     return nFaceSets;
 }
 
 
-// Remove points not used by any face or points used by only two faces where
-// the edges are in line
-Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::mergeEdges
+// Remove points. pointCanBeDeleted is parallel synchronised.
+Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::doRemovePoints
 (
-    const scalar minCos
+    removePoints& pointRemover,
+    const boolList& pointCanBeDeleted
 )
 {
-    // Point removal analysis engine
-    removePoints pointRemover(mesh_);
+    // Topology changes container
+    polyTopoChange meshMod(mesh_);
 
-    // Count usage of points
-    boolList pointCanBeDeleted;
-    label nRemove = pointRemover.countPointUsage(minCos, pointCanBeDeleted);
+    pointRemover.setRefinement(pointCanBeDeleted, meshMod);
 
-    Info<< "Removing " << nRemove
-        << " straight edge points." << endl;
+    // Change the mesh (no inflation)
+    autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
 
-    autoPtr<mapPolyMesh> map;
+    // Update fields
+    mesh_.updateMesh(map);
 
-    if (nRemove > 0)
+    // Move mesh (since morphing does not do this)
+    if (map().hasMotionPoints())
     {
-        // Save my local faces that will change. These changed faces might
-        // cause a shift in the cell centre which needs to be retested.
-        // Have to do this before changing mesh since point will be removed.
-        labelHashSet retestOldFaces(nRemove / Pstream::nProcs());
-
-        {
-            const faceList& faces = mesh_.faces();
-
-            forAll(faces, faceI)
-            {
-                const face& f = faces[faceI];
-
-                forAll(f, fp)
-                {
-                    if (pointCanBeDeleted[f[fp]])
-                    {
-                        retestOldFaces.insert(faceI);
-                        break;
-                    }
-                }
-            }
-        }
-
-        // Topology changes container
-        polyTopoChange meshMod(mesh_);
-
-        pointRemover.setRefinement(pointCanBeDeleted, meshMod);
-
-        // Change the mesh (no inflation)
-        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();
-        }
-
-        pointRemover.updateMesh(map);
-
-        // Get the kept faces that need to be recalculated.
-        labelHashSet retestFaces(6*retestOldFaces.size());
-
-        const cellList& cells = mesh_.cells();
-
-        forAllConstIter(labelHashSet, retestOldFaces, iter)
-        {
-            label faceI = map().reverseFaceMap()[iter.key()];
-
-            const cell& ownFaces = cells[mesh_.faceOwner()[faceI]];
-
-            forAll(ownFaces, i)
-            {
-                retestFaces.insert(ownFaces[i]);
-            }
-
-            if (mesh_.isInternalFace(faceI))
-            {
-                const cell& neiFaces = cells[mesh_.faceNeighbour()[faceI]];
-
-                forAll(neiFaces, i)
-                {
-                    retestFaces.insert(neiFaces[i]);
-                }
-            }
-        }
-        updateMesh(map, retestFaces.toc());
+        mesh_.movePoints(map().preMotionPoints());
     }
+    else
+    {
+        // Delete mesh volumes.
+        mesh_.clearOut();
+    }
+
+    // Reset the instance for if in overwrite mode
+    mesh_.setInstance(timeName());
+
+    pointRemover.updateMesh(map);
+    updateMesh(map, labelList(0));
 
     return map;
 }
 
 
+// Restore faces (which contain removed points)
+Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::doRestorePoints
+(
+    removePoints& pointRemover,
+    const labelList& facesToRestore
+)
+{
+    // Topology changes container
+    polyTopoChange meshMod(mesh_);
+
+    // Determine sets of points and faces to restore
+    labelList localFaces, localPoints;
+    pointRemover.getUnrefimentSet
+    (
+        facesToRestore,
+        localFaces,
+        localPoints
+    );
+
+    // Undo the changes on the faces that are in error.
+    pointRemover.setUnrefinement
+    (
+        localFaces,
+        localPoints,
+        meshMod
+    );
+
+    // Change the mesh (no inflation)
+    autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
+
+    // Update fields
+    mesh_.updateMesh(map);
+
+    // Move mesh (since morphing does not do this)
+    if (map().hasMotionPoints())
+    {
+        mesh_.movePoints(map().preMotionPoints());
+    }
+    else
+    {
+        // Delete mesh volumes.
+        mesh_.clearOut();
+    }
+
+    // Reset the instance for if in overwrite mode
+    mesh_.setInstance(timeName());
+
+    pointRemover.updateMesh(map);
+    updateMesh(map, labelList(0));
+
+    return map;
+}
+
+
+// Collect all faces that are both in candidateFaces and in the set.
+// If coupled face also collects the coupled face.
+Foam::labelList Foam::meshRefinement::collectFaces
+(
+    const labelList& candidateFaces,
+    const labelHashSet& set
+) const
+{
+    // Has face been selected?
+    boolList selected(mesh_.nFaces(), false);
+
+    forAll(candidateFaces, i)
+    {
+        label faceI = candidateFaces[i];
+
+        if (set.found(faceI))
+        {
+            selected[faceI] = true;
+        }
+    }
+    syncTools::syncFaceList
+    (
+        mesh_,
+        selected,
+        orEqOp<bool>()      // combine operator
+    );
+
+    labelList selectedFaces(findIndices(selected, true));
+
+    return selectedFaces;
+}
+
+
+// Pick up faces of cells of faces in set.
+Foam::labelList Foam::meshRefinement::growFaceCellFace
+(
+    const labelHashSet& set
+) const
+{
+    boolList selected(mesh_.nFaces(), false);
+
+    forAllConstIter(faceSet, set, iter)
+    {
+        label faceI = iter.key();
+
+        label own = mesh_.faceOwner()[faceI];
+
+        const cell& ownFaces = mesh_.cells()[own];
+        forAll(ownFaces, ownFaceI)
+        {
+            selected[ownFaces[ownFaceI]] = true;
+        }
+
+        if (mesh_.isInternalFace(faceI))
+        {
+            label nbr = mesh_.faceNeighbour()[faceI];
+
+            const cell& nbrFaces = mesh_.cells()[nbr];
+            forAll(nbrFaces, nbrFaceI)
+            {
+                selected[nbrFaces[nbrFaceI]] = true;
+            }
+        }
+    }
+    syncTools::syncFaceList
+    (
+        mesh_,
+        selected,
+        orEqOp<bool>()      // combine operator
+    );
+    return findIndices(selected, true);
+}
+
+
+// Remove points not used by any face or points used by only two faces where
+// the edges are in line
+Foam::label Foam::meshRefinement::mergeEdgesUndo
+(
+    const scalar minCos,
+    const dictionary& motionDict
+)
+{
+    Info<< nl
+        << "Merging all points on surface that" << nl
+        << "- are used by only two boundary faces and" << nl
+        << "- make an angle with a cosine of more than " << minCos
+        << "." << nl << endl;
+
+    // Point removal analysis engine with undo
+    removePoints pointRemover(mesh_, true);
+
+    // Count usage of points
+    boolList pointCanBeDeleted;
+    label nRemove = pointRemover.countPointUsage(minCos, pointCanBeDeleted);
+
+    if (nRemove > 0)
+    {
+        Info<< "Removing " << nRemove
+            << " straight edge points ..." << nl << endl;
+
+        // Remove points
+        // ~~~~~~~~~~~~~
+
+        doRemovePoints(pointRemover, pointCanBeDeleted);
+
+
+        for (label iteration = 0; iteration < 100; iteration++)
+        {
+            Info<< nl
+                << "Undo iteration " << iteration << nl
+                << "----------------" << endl;
+
+
+            // Check mesh for errors
+            // ~~~~~~~~~~~~~~~~~~~~~
+
+            faceSet errorFaces
+            (
+                mesh_,
+                "errorFaces",
+                mesh_.nFaces()-mesh_.nInternalFaces()
+            );
+            bool hasErrors = motionSmoother::checkMesh
+            (
+                false,  // report
+                mesh_,
+                motionDict,
+                errorFaces
+            );
+            //if (checkEdgeConnectivity)
+            //{
+            //    Info<< "Checking edge-face connectivity (duplicate faces"
+            //        << " or non-consecutive shared vertices)" << endl;
+            //
+            //    label nOldSize = errorFaces.size();
+            //
+            //    hasErrors =
+            //        mesh_.checkFaceFaces
+            //        (
+            //            false,
+            //            &errorFaces
+            //        )
+            //     || hasErrors;
+            //
+            //    Info<< "Detected additional "
+            //        << returnReduce(errorFaces.size()-nOldSize,sumOp<label>())
+            //        << " faces with illegal face-face connectivity"
+            //        << endl;
+            //}
+
+            if (!hasErrors)
+            {
+                break;
+            }
+
+            if (debug)
+            {
+                errorFaces.instance() = timeName();
+                Pout<< "**Writing all faces in error to faceSet "
+                    << errorFaces.objectPath() << nl << endl;
+                errorFaces.write();
+            }
+
+            labelList masterErrorFaces
+            (
+                collectFaces
+                (
+                    pointRemover.savedFaceLabels(),
+                    errorFaces
+                )
+            );
+
+            label n = returnReduce(masterErrorFaces.size(), sumOp<label>());
+
+            Info<< "Detected " << n
+                << " error faces on boundaries that have been merged."
+                << " These will be restored to their original faces." << nl
+                << endl;
+
+            if (n == 0)
+            {
+                if (hasErrors)
+                {
+                    Info<< "Detected "
+                        << returnReduce(errorFaces.size(), sumOp<label>())
+                        << " error faces in mesh."
+                        << " Restoring neighbours of faces in error." << nl
+                        << endl;
+
+                    labelList expandedErrorFaces
+                    (
+                        growFaceCellFace
+                        (
+                            errorFaces
+                        )
+                    );
+
+                    doRestorePoints(pointRemover, expandedErrorFaces);
+                }
+
+                break;
+            }
+
+            doRestorePoints(pointRemover, masterErrorFaces);
+        }
+
+        if (debug)
+        {
+            const_cast<Time&>(mesh_.time())++;
+            Pout<< "Writing merged-edges mesh to time "
+                << timeName() << nl << endl;
+            write();
+        }
+    }
+    else
+    {
+        Info<< "No straight edges simplified and no points removed ..." << endl;
+    }
+
+    return nRemove;
+}
+
+
 // ************************************************************************* //

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -202,7 +202,7 @@ Foam::Map<Foam::label> Foam::meshRefinement::findEdgeConnectedProblemCells
     if (debug)
     {
         faceSet fSet(mesh_, "edgeConnectedFaces", candidateFaces);
-        fSet.instance() = mesh_.time().timeName();
+        fSet.instance() = timeName();
         Pout<< "Writing " << fSet.size()
             << " with problematic topology to faceSet "
             << fSet.objectPath() << endl;
@@ -265,7 +265,7 @@ Foam::Map<Foam::label> Foam::meshRefinement::findEdgeConnectedProblemCells
 
     if (debug)
     {
-        perpFaces.instance() = mesh_.time().timeName();
+        perpFaces.instance() = timeName();
         Pout<< "Writing " << perpFaces.size()
             << " faces that are perpendicular to the surface to set "
             << perpFaces.objectPath() << endl;
@@ -427,6 +427,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
 
     // Count of faces marked for baffling
     label nBaffleFaces = 0;
+    PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
 
     // Count of faces not baffled since would not cause a collapse
     label nPrevented = 0;
@@ -483,7 +484,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
         if (debug)
         {
             cellSet problemCellSet(mesh_, "problemCells", problemCells.toc());
-            problemCellSet.instance() = mesh_.time().timeName();
+            problemCellSet.instance() = timeName();
             Pout<< "Writing " << problemCellSet.size()
                 << " cells that are edge connected to coarser cell to set "
                 << problemCellSet.objectPath() << endl;
@@ -918,7 +919,10 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
                         else
                         {
                             facePatch[faceI] = getBafflePatch(facePatch, faceI);
-                            nBaffleFaces++;
+                            if (isMasterFace[faceI])
+                            {
+                                nBaffleFaces++;
+                            }
 
                             // Do NOT update boundary data since this would grow
                             // blocked faces across gaps.

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

@@ -28,6 +28,7 @@ License
 #include "syncTools.H"
 #include "Time.H"
 #include "refinementSurfaces.H"
+#include "refinementFeatures.H"
 #include "shellSurfaces.H"
 #include "faceSet.H"
 #include "decompositionMethod.H"
@@ -53,6 +54,11 @@ Foam::labelList Foam::meshRefinement::getChangedFaces
     const polyMesh& mesh = map.mesh();
 
     labelList changedFaces;
+
+    // For reporting: number of masterFaces changed
+    label nMasterChanged = 0;
+    PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh));
+
     {
         // Mark any face on a cell which has been added or changed
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -182,36 +188,26 @@ Foam::labelList Foam::meshRefinement::getChangedFaces
         // Now we have in changedFace marked all affected faces. Pack.
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-        label nChanged = 0;
+        changedFaces = findIndices(changedFace, true);
+
+        // Count changed master faces.
+        nMasterChanged = 0;
 
         forAll(changedFace, faceI)
         {
-            if (changedFace[faceI])
+            if (changedFace[faceI] && isMasterFace[faceI])
             {
-                nChanged++;
+                nMasterChanged++;
             }
         }
 
-        changedFaces.setSize(nChanged);
-        nChanged = 0;
-
-        forAll(changedFace, faceI)
-        {
-            if (changedFace[faceI])
-            {
-                changedFaces[nChanged++] = faceI;
-            }
-        }
     }
 
-    label nChangedFaces = changedFaces.size();
-    reduce(nChangedFaces, sumOp<label>());
-
     if (debug)
     {
         Pout<< "getChangedFaces : Detected "
             << " local:" << changedFaces.size()
-            << " global:" << nChangedFaces
+            << " global:" << returnReduce(nMasterChanged, sumOp<label>())
             << " changed faces out of " << mesh.globalData().nTotalFaces()
             << endl;
 
@@ -252,8 +248,6 @@ bool Foam::meshRefinement::markForRefine
 Foam::label Foam::meshRefinement::markFeatureRefinement
 (
     const point& keepPoint,
-    const PtrList<featureEdgeMesh>& featureMeshes,
-    const labelList& featureLevels,
     const label nAllowRefine,
 
     labelList& refineCell,
@@ -287,9 +281,11 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
 
     if (cellI != -1)
     {
-        forAll(featureMeshes, featI)
+        forAll(features_, featI)
         {
-            const featureEdgeMesh& featureMesh = featureMeshes[featI];
+            const featureEdgeMesh& featureMesh = features_[featI];
+            const label featureLevel = features_.levels()[featI];
+
             const labelListList& pointEdges = featureMesh.pointEdges();
 
             forAll(pointEdges, pointI)
@@ -315,7 +311,7 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
                             tetFaceI,
                             tetPtI,
                             featureMesh.points()[pointI],   // endpos
-                            featureLevels[featI],           // level
+                            featureLevel,                   // level
                             featI,                          // featureMesh
                             pointI                          // end point
                         )
@@ -339,11 +335,11 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
     maxFeatureLevel = -1;
 
     // Whether edge has been visited.
-    List<PackedBoolList> featureEdgeVisited(featureMeshes.size());
+    List<PackedBoolList> featureEdgeVisited(features_.size());
 
-    forAll(featureMeshes, featI)
+    forAll(features_, featI)
     {
-        featureEdgeVisited[featI].setSize(featureMeshes[featI].edges().size());
+        featureEdgeVisited[featI].setSize(features_[featI].edges().size());
         featureEdgeVisited[featI] = 0u;
     }
 
@@ -359,7 +355,7 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
             label featI = tp.i();
             label pointI = tp.j();
 
-            const featureEdgeMesh& featureMesh = featureMeshes[featI];
+            const featureEdgeMesh& featureMesh = features_[featI];
             const labelList& pEdges = featureMesh.pointEdges()[pointI];
 
             // Particle now at pointI. Check connected edges to see which one
@@ -1068,9 +1064,6 @@ Foam::labelList Foam::meshRefinement::refineCandidates
     const point& keepPoint,
     const scalar curvature,
 
-    const PtrList<featureEdgeMesh>& featureMeshes,
-    const labelList& featureLevels,
-
     const bool featureRefinement,
     const bool internalRefinement,
     const bool surfaceRefinement,
@@ -1134,8 +1127,6 @@ Foam::labelList Foam::meshRefinement::refineCandidates
             label nFeatures = markFeatureRefinement
             (
                 keepPoint,
-                featureMeshes,
-                featureLevels,
                 nAllowRefine,
 
                 refineCell,
@@ -1251,6 +1242,9 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::refine
         mesh_.clearOut();
     }
 
+    // Reset the instance for if in overwrite mode
+    mesh_.setInstance(timeName());
+
     // Update intersection info
     updateMesh(map, getChangedFaces(map, cellsToRefine));
 

src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C

@@ -0,0 +1,233 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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 "refinementFeatures.H"
+#include "Time.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::refinementFeatures::refinementFeatures
+(
+    const objectRegistry& io,
+    const PtrList<dictionary>& featDicts
+)
+:
+    PtrList<featureEdgeMesh>(featDicts.size()),
+    levels_(featDicts.size()),
+    edgeTrees_(featDicts.size()),
+    pointTrees_(featDicts.size())
+{
+    // Read features
+
+    forAll(featDicts, i)
+    {
+        const dictionary& dict = featDicts[i];
+
+        fileName featFileName(dict.lookup("file"));
+
+        set
+        (
+            i,
+            new featureEdgeMesh
+            (
+                IOobject
+                (
+                    featFileName,                       // name
+                    io.time().constant(),               // instance
+                    "triSurface",                       // local
+                    io.time(),                          // registry
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE,
+                    false
+                )
+            )
+        );
+
+        const featureEdgeMesh& eMesh = operator[](i);
+
+        //eMesh.mergePoints(meshRefiner_.mergeDistance());
+        levels_[i] = readLabel(dict.lookup("level"));
+
+        Info<< "Refinement level " << levels_[i]
+            << " for all cells crossed by feature " << featFileName
+            << " (" << eMesh.points().size() << " points, "
+            << eMesh.edges().size() << " edges)." << endl;
+    }
+
+
+    // Search engines
+
+    forAll(*this, i)
+    {
+        const featureEdgeMesh& eMesh = operator[](i);
+        const pointField& points = eMesh.points();
+        const edgeList& edges = eMesh.edges();
+
+        // Calculate bb of all points
+        const treeBoundBox bb(points);
+
+        edgeTrees_.set
+        (
+            i,
+            new indexedOctree<treeDataEdge>
+            (
+                treeDataEdge
+                (
+                    false,                  // do not cache bb
+                    edges,
+                    points,
+                    identity(edges.size())
+                ),
+                bb,     // overall search domain
+                8,      // maxLevel
+                10,     // leafsize
+                3.0     // duplicity
+            )
+        );
+
+
+        // Detect feature points from edges.
+        const labelListList& pointEdges = eMesh.pointEdges();
+        DynamicList<label> featurePoints;
+        forAll(pointEdges, pointI)
+        {
+            if (pointEdges[pointI].size() > 2)
+            {
+                featurePoints.append(pointI);
+            }
+        }
+
+        Info<< "Detected " << featurePoints.size()
+            << " featurePoints out of " << points.size() << endl;
+
+        pointTrees_.set
+        (
+            i,
+            new indexedOctree<treeDataPoint>
+            (
+                treeDataPoint(points, featurePoints),
+                bb,     // overall search domain
+                8,      // maxLevel
+                10,     // leafsize
+                3.0     // duplicity
+            )
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::refinementFeatures::findNearestEdge
+(
+    const pointField& samples,
+    const scalarField& nearestDistSqr,
+    labelList& nearFeature,
+    List<pointIndexHit>& nearInfo
+) const
+{
+    nearFeature.setSize(samples.size());
+    nearFeature = -1;
+    nearInfo.setSize(samples.size());
+
+    forAll(edgeTrees_, featI)
+    {
+        const indexedOctree<treeDataEdge>& tree = edgeTrees_[featI];
+        forAll(samples, sampleI)
+        {
+            const point& sample = samples[sampleI];
+
+            scalar distSqr;
+            if (nearInfo[sampleI].hit())
+            {
+                distSqr = magSqr(nearInfo[sampleI].hitPoint()-sample);
+            }
+            else
+            {
+                distSqr = nearestDistSqr[sampleI];
+            }
+
+            pointIndexHit info = tree.findNearest(sample, distSqr);
+
+            if (info.hit())
+            {
+                nearInfo[sampleI] = info;
+                nearFeature[sampleI] = featI;
+            }
+        }
+    }
+}
+
+
+void Foam::refinementFeatures::findNearestPoint
+(
+    const pointField& samples,
+    const scalarField& nearestDistSqr,
+    labelList& nearFeature,
+    labelList& nearIndex
+) const
+{
+    nearFeature.setSize(samples.size());
+    nearFeature = -1;
+    nearIndex.setSize(samples.size());
+    nearIndex = -1;
+
+    forAll(pointTrees_, featI)
+    {
+        const indexedOctree<treeDataPoint>& tree = pointTrees_[featI];
+        forAll(samples, sampleI)
+        {
+            const point& sample = samples[sampleI];
+
+            scalar distSqr;
+            if (nearFeature[sampleI] != -1)
+            {
+                label nearFeatI = nearFeature[sampleI];
+                const indexedOctree<treeDataPoint>& nearTree =
+                    pointTrees_[nearFeatI];
+                label featPointI =
+                    nearTree.shapes().pointLabels()[nearIndex[sampleI]];
+                const point& featPt =
+                    operator[](nearFeatI).points()[featPointI];
+                distSqr = magSqr(featPt-sample);
+            }
+            else
+            {
+                distSqr = nearestDistSqr[sampleI];
+            }
+
+            pointIndexHit info = tree.findNearest(sample, distSqr);
+
+            if (info.hit())
+            {
+                nearFeature[sampleI] = featI;
+                nearIndex[sampleI] = info.index();
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H

@@ -0,0 +1,139 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 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::refinementFeatures
+
+Description
+    Encapsulates queries for features.
+
+SourceFiles
+    refinementFeatures.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef refinementFeatures_H
+#define refinementFeatures_H
+
+#include "featureEdgeMesh.H"
+#include "indexedOctree.H"
+#include "treeDataEdge.H"
+#include "treeDataPoint.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class refinementFeatures Declaration
+\*---------------------------------------------------------------------------*/
+
+class refinementFeatures
+:
+    public PtrList<featureEdgeMesh>
+{
+private:
+
+    // Private data
+
+        //- Refinement levels
+        labelList levels_;
+
+        //- Edge
+        PtrList<indexedOctree<treeDataEdge> > edgeTrees_;
+
+        //- Features points
+        PtrList<indexedOctree<treeDataPoint> > pointTrees_;
+
+
+    // Private Member Functions
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        refinementFeatures
+        (
+            const objectRegistry& io,
+            const PtrList<dictionary>& featDicts
+        );
+
+
+    // Member Functions
+
+        // Access
+
+            const labelList& levels() const
+            {
+                return levels_;
+            }
+
+            const PtrList<indexedOctree<treeDataEdge> >& edgeTrees() const
+            {
+                return edgeTrees_;
+            }
+
+            const PtrList<indexedOctree<treeDataPoint> >& pointTrees() const
+            {
+                return pointTrees_;
+            }
+
+
+        // Query
+
+            //- Find nearest point on nearest feature edge
+            void findNearestEdge
+            (
+                const pointField& samples,
+                const scalarField& nearestDistSqr,
+                labelList& nearFeature,
+                List<pointIndexHit>& nearInfo
+            ) const;
+
+            //- Find nearest feature point. Is an index into feature points
+            //  which itself is an index into the edgeMesh points.
+            //  So the point index is
+            //      pointTrees()[nearFeature].shapes().pointLabels()[nearIndex]
+            //  Wip.
+            void findNearestPoint
+            (
+                const pointField& samples,
+                const scalarField& nearestDistSqr,
+                labelList& nearFeature,
+                labelList& nearIndex
+            ) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/sets/cellSources/nearestToCell/nearestToCell.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -55,12 +55,9 @@ Foam::topoSetSource::addToUsageTable Foam::nearestToCell::usage_
 
 void Foam::nearestToCell::combine(topoSet& set, const bool add) const
 {
-    // Construct search engine withouth tet decomposition.
-    meshSearch queryMesh(mesh_, false);
-
     forAll(points_, pointI)
     {
-        addOrDelete(set, queryMesh.findNearestCell(points_[pointI]), add);
+        addOrDelete(set, mesh_.findNearestCell(points_[pointI]), add);
     }
 }
 

src/meshTools/sets/topoSets/faceSet.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -140,7 +140,7 @@ void faceSet::sync(const polyMesh& mesh)
     }
 
     reduce(nAdded, sumOp<label>());
-    if (nAdded > 0)
+    if (debug && nAdded > 0)
     {
         Info<< "Added an additional " << nAdded
             << " faces on coupled patches. "

src/regionModels/regionModel/regionModel1D/regionModel1D.C

@@ -114,6 +114,8 @@ void Foam::regionModels::regionModel1D::initialise()
             boundaryFaceCells_[localPyrolysisFaceI].transfer(cellIDs);
 
             localPyrolysisFaceI++;
+
+            nLayers_ = nCells;
         }
     }
 
@@ -266,6 +268,7 @@ Foam::regionModels::regionModel1D::regionModel1D(const fvMesh& mesh)
     boundaryFaceFaces_(),
     boundaryFaceCells_(),
     boundaryFaceOppositeFace_(),
+    nLayers_(0),
     nMagSfPtr_(NULL),
     moveMesh_(false)
 {}
@@ -283,6 +286,7 @@ Foam::regionModels::regionModel1D::regionModel1D
     boundaryFaceFaces_(regionMesh().nCells()),
     boundaryFaceCells_(regionMesh().nCells()),
     boundaryFaceOppositeFace_(regionMesh().nCells()),
+    nLayers_(0),
     nMagSfPtr_(NULL),
     moveMesh_(true)
 {

src/regionModels/regionModel/regionModel1D/regionModel1D.H

@@ -88,6 +88,9 @@ protected:
             //- Global boundary face IDs oppossite coupled patch
             labelList boundaryFaceOppositeFace_;
 
+            //- Number of layers in the region
+            label nLayers_;
+
 
         // Geometry
 
@@ -115,7 +118,7 @@ protected:
 public:
 
     //- Runtime type information
-    TypeName("regionModel");
+    TypeName("regionModel1D");
 
 
     // Constructors
@@ -152,6 +155,9 @@ public:
                 //- Return the global boundary face IDs oppossite coupled patch
                 inline const labelList& boundaryFaceOppositeFace() const;
 
+                //- Return the number of layers in the region
+                inline label nLayers() const;
+
 
             // Geometry
 

src/regionModels/regionModel/regionModel1D/regionModel1DI.H

@@ -49,6 +49,12 @@ Foam::regionModels::regionModel1D::boundaryFaceOppositeFace() const
 }
 
 
+inline Foam::label Foam::regionModels::regionModel1D::nLayers() const
+{
+    return nLayers_;
+}
+
+
 inline const Foam::surfaceScalarField&
 Foam::regionModels::regionModel1D::nMagSf() const
 {

src/sampling/sampledSurface/sampledPatch/sampledPatch.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -46,15 +46,14 @@ Foam::sampledPatch::sampledPatch
 (
     const word& name,
     const polyMesh& mesh,
-    const word& patchName,
+    const wordReList& patchNames,
     const bool triangulate
 )
 :
     sampledSurface(name, mesh),
-    patchName_(patchName),
+    patchNames_(patchNames),
     triangulate_(triangulate),
-    needsUpdate_(true),
-    patchFaceLabels_(0)
+    needsUpdate_(true)
 {}
 
 
@@ -66,10 +65,9 @@ Foam::sampledPatch::sampledPatch
 )
 :
     sampledSurface(name, mesh, dict),
-    patchName_(dict.lookup("patchName")),
+    patchNames_(dict.lookup("patches")),
     triangulate_(dict.lookupOrDefault("triangulate", false)),
-    needsUpdate_(true),
-    patchFaceLabels_(0)
+    needsUpdate_(true)
 {}
 
 
@@ -81,6 +79,20 @@ Foam::sampledPatch::~sampledPatch()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+const Foam::labelList& Foam::sampledPatch::patchIDs() const
+{
+    if (patchIDs_.empty())
+    {
+        patchIDs_ = mesh().boundaryMesh().patchSet
+        (
+            patchNames_,
+            false
+        ).sortedToc();
+    }
+    return patchIDs_;
+}
+
+
 bool Foam::sampledPatch::needsUpdate() const
 {
     return needsUpdate_;
@@ -97,7 +109,10 @@ bool Foam::sampledPatch::expire()
 
     sampledSurface::clearGeom();
     MeshStorage::clear();
+    patchIDs_.clear();
+    patchIndex_.clear();
     patchFaceLabels_.clear();
+    patchStart_.clear();
 
     needsUpdate_ = true;
     return true;
@@ -111,31 +126,67 @@ bool Foam::sampledPatch::update()
         return false;
     }
 
-    const label patchI = mesh().boundaryMesh().findPatchID(patchName_);
-
-    if (patchI != -1)
+    label sz = 0;
+    forAll(patchIDs(), i)
     {
-        const polyPatch& p = mesh().boundaryMesh()[patchI];
-        this->storedPoints() = p.localPoints();
-        this->storedFaces()  = p.localFaces();
+        label patchI = patchIDs()[i];
+        const polyPatch& pp = mesh().boundaryMesh()[patchI];
 
-        // an identity map
-        patchFaceLabels_.setSize(faces().size());
-        forAll(patchFaceLabels_, i)
+        if (isA<emptyPolyPatch>(pp))
         {
-            patchFaceLabels_[i] = i;
+            FatalErrorIn("sampledPatch::update()")
+                << "Cannot sample an empty patch. Patch " << pp.name()
+                << exit(FatalError);
         }
 
-        // triangulate uses remapFaces()
-        // - this is somewhat less efficient since it recopies the faces
-        // that we just created, but we probably don't want to do this
-        // too often anyhow.
-        if (triangulate_)
+        sz += pp.size();
+    }
+
+    // For every face (or triangle) the originating patch and local face in the
+    // patch.
+    patchIndex_.setSize(sz);
+    patchFaceLabels_.setSize(sz);
+    patchStart_.setSize(patchIDs().size());
+    labelList meshFaceLabels(sz);
+
+    sz = 0;
+
+    forAll(patchIDs(), i)
+    {
+        label patchI = patchIDs()[i];
+
+        patchStart_[i] = sz;
+
+        const polyPatch& pp = mesh().boundaryMesh()[patchI];
+
+        forAll(pp, j)
         {
-            MeshStorage::triangulate();
+            patchIndex_[sz] = i;
+            patchFaceLabels_[sz] = j;
+            meshFaceLabels[sz] = pp.start()+j;
+            sz++;
         }
     }
 
+    indirectPrimitivePatch allPatches
+    (
+        IndirectList<face>(mesh().faces(), meshFaceLabels),
+        mesh().points()
+    );
+
+    this->storedPoints() = allPatches.localPoints();
+    this->storedFaces()  = allPatches.localFaces();
+
+
+    // triangulate uses remapFaces()
+    // - this is somewhat less efficient since it recopies the faces
+    // that we just created, but we probably don't want to do this
+    // too often anyhow.
+    if (triangulate_)
+    {
+        MeshStorage::triangulate();
+    }
+
     if (debug)
     {
         print(Pout);
@@ -148,10 +199,7 @@ bool Foam::sampledPatch::update()
 
 
 // remap action on triangulation
-void Foam::sampledPatch::remapFaces
-(
-    const labelUList& faceMap
-)
+void Foam::sampledPatch::remapFaces(const labelUList& faceMap)
 {
     // recalculate the cells cut
     if (&faceMap && faceMap.size())
@@ -161,11 +209,27 @@ void Foam::sampledPatch::remapFaces
         (
             UIndirectList<label>(patchFaceLabels_, faceMap)
         );
+        patchIndex_ = labelList
+        (
+            UIndirectList<label>(patchIndex_, faceMap)
+        );
+
+        // Redo patchStart.
+        if (patchIndex_.size() > 0)
+        {
+            patchStart_[patchIndex_[0]] = 0;
+            for (label i = 1; i < patchIndex_.size(); i++)
+            {
+                if (patchIndex_[i] != patchIndex_[i-1])
+                {
+                    patchStart_[patchIndex_[i]] = i;
+                }
+            }
+        }
     }
 }
 
 
-
 Foam::tmp<Foam::scalarField> Foam::sampledPatch::sample
 (
     const volScalarField& vField
@@ -257,7 +321,7 @@ Foam::tmp<Foam::tensorField> Foam::sampledPatch::interpolate
 void Foam::sampledPatch::print(Ostream& os) const
 {
     os  << "sampledPatch: " << name() << " :"
-        << "  patch:" << patchName()
+        << "  patches:" << patchNames()
         << "  faces:" << faces().size()
         << "  points:" << points().size();
 }

src/sampling/sampledSurface/sampledPatch/sampledPatch.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,7 +25,7 @@ Class
     Foam::sampledPatch
 
 Description
-    A sampledSurface on a patch. Non-triangulated by default.
+    A sampledSurface on patches. Non-triangulated by default.
 
 SourceFiles
     sampledPatch.C
@@ -57,8 +57,11 @@ class sampledPatch
 
     // Private data
 
-        //- Name of patch
-        const word patchName_;
+        //- Name of patches
+        const wordReList patchNames_;
+
+        //- Corresponding patchIDs
+        mutable labelList patchIDs_;
 
         //- Triangulated faces or keep faces as is
         bool triangulate_;
@@ -66,11 +69,16 @@ class sampledPatch
         //- Track if the surface needs an update
         mutable bool needsUpdate_;
 
-        //- Local patch face labels
+        //- For every face (or triangle) the originating patch
+        labelList patchIndex_;
+
+        //- For every face (or triangle) the index in the originating patch
         labelList patchFaceLabels_;
 
-    // Private Member Functions
+        //- Start indices (in patchFaceLabels_) of patches
+        labelList patchStart_;
 
+    // Private Member Functions
 
         //- sample field on faces
         template <class Type>
@@ -79,15 +87,34 @@ class sampledPatch
             const GeometricField<Type, fvPatchField, volMesh>& vField
         ) const;
 
-
         template <class Type>
         tmp<Field<Type> >
         interpolateField(const interpolation<Type>&) const;
 
-
         //- remap action on triangulation or cleanup
         virtual void remapFaces(const labelUList& faceMap);
 
+
+protected:
+
+        const wordReList& patchNames() const
+        {
+            return patchNames_;
+        }
+
+        const labelList& patchIDs() const;
+
+        const labelList& patchStart() const
+        {
+            return patchStart_;
+        }
+
+        const labelList& patchFaceLabels() const
+        {
+            return patchFaceLabels_;
+        }
+
+
 public:
 
     //- Runtime type information
@@ -101,7 +128,7 @@ public:
         (
             const word& name,
             const polyMesh& mesh,
-            const word& patchName,
+            const wordReList& patchNames,
             const bool triangulate = false
         );
 
@@ -133,21 +160,6 @@ public:
         virtual bool update();
 
 
-        const word patchName() const
-        {
-            return patchName_;
-        }
-
-        label patchIndex() const
-        {
-            return mesh().boundaryMesh().findPatchID(patchName_);
-        }
-
-        const labelList& patchFaceLabels() const
-        {
-            return patchFaceLabels_;
-        }
-
         //- Points of surface
         virtual const pointField& points() const
         {

src/sampling/sampledSurface/sampledPatch/sampledPatchTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,15 +37,11 @@ Foam::sampledPatch::sampleField
     // One value per face
     tmp<Field<Type> > tvalues(new Field<Type>(patchFaceLabels_.size()));
     Field<Type>& values = tvalues();
-
-    if (patchIndex() != -1)
+    forAll(patchFaceLabels_, i)
     {
-        const Field<Type>& bField = vField.boundaryField()[patchIndex()];
-
-        forAll(patchFaceLabels_, elemI)
-        {
-            values[elemI] = bField[patchFaceLabels_[elemI]];
-        }
+        label patchI = patchIDs_[patchIndex_[i]];
+        const Field<Type>& bField = vField.boundaryField()[patchI];
+        values[i] = bField[patchFaceLabels_[i]];
     }
 
     return tvalues;
@@ -63,34 +59,33 @@ Foam::sampledPatch::interpolateField
     tmp<Field<Type> > tvalues(new Field<Type>(points().size()));
     Field<Type>& values = tvalues();
 
-    if (patchIndex() != -1)
+    const labelList& own = mesh().faceOwner();
+
+    boolList pointDone(points().size(), false);
+
+    forAll(faces(), cutFaceI)
     {
-        const polyPatch& patch = mesh().boundaryMesh()[patchIndex()];
-        const labelList& own = mesh().faceOwner();
+        label patchI = patchIDs_[patchIndex_[cutFaceI]];
+        const polyPatch& pp = mesh().boundaryMesh()[patchI];
+        label patchFaceI = patchFaceLabels()[cutFaceI];
+        const face& f = faces()[cutFaceI];
 
-        boolList pointDone(points().size(), false);
-
-        forAll(faces(), cutFaceI)
+        forAll(f, faceVertI)
         {
-            const face& f = faces()[cutFaceI];
+            label pointI = f[faceVertI];
 
-            forAll(f, faceVertI)
+            if (!pointDone[pointI])
             {
-                label pointI = f[faceVertI];
+                label faceI = patchFaceI + pp.start();
+                label cellI = own[faceI];
 
-                if (!pointDone[pointI])
-                {
-                    label faceI = patchFaceLabels()[cutFaceI] + patch.start();
-                    label cellI = own[faceI];
-
-                    values[pointI] = interpolator.interpolate
-                    (
-                        points()[pointI],
-                        cellI,
-                        faceI
-                    );
-                    pointDone[pointI] = true;
-                }
+                values[pointI] = interpolator.interpolate
+                (
+                    points()[pointI],
+                    cellI,
+                    faceI
+                );
+                pointDone[pointI] = true;
             }
         }
     }

src/sampling/sampledSurface/sampledPatchInternalField/sampledPatchInternalField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,15 +56,27 @@ Foam::sampledPatchInternalField::sampledPatchInternalField
 )
 :
     sampledPatch(name, mesh, dict),
-    directMappedPatchBase
-    (
-        mesh.boundaryMesh()[sampledPatch::patchIndex()],
-        mesh.name(),                        // sampleRegion
-        directMappedPatchBase::NEARESTCELL, // sampleMode
-        word::null,                         // samplePatch
-        -readScalar(dict.lookup("distance"))
-    )
-{}
+    mappers_(patchIDs().size())
+{
+    const scalar distance = readScalar(dict.lookup("distance"));
+
+    forAll(patchIDs(), i)
+    {
+        label patchI = patchIDs()[i];
+        mappers_.set
+        (
+            i,
+            new directMappedPatchBase
+            (
+                mesh.boundaryMesh()[patchI],
+                mesh.name(),                        // sampleRegion
+                directMappedPatchBase::NEARESTCELL, // sampleMode
+                word::null,                         // samplePatch
+                -distance                           // sample inside my domain
+            )
+        );
+    }
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -167,7 +179,7 @@ Foam::tmp<Foam::tensorField> Foam::sampledPatchInternalField::interpolate
 void Foam::sampledPatchInternalField::print(Ostream& os) const
 {
     os  << "sampledPatchInternalField: " << name() << " :"
-        << "  patch:" << patchName()
+        << "  patches:" << patchNames()
         << "  faces:" << faces().size()
         << "  points:" << points().size();
 }

src/sampling/sampledSurface/sampledPatchInternalField/sampledPatchInternalField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -54,9 +54,13 @@ namespace Foam
 
 class sampledPatchInternalField
 :
-    public sampledPatch,
-    public directMappedPatchBase
+    public sampledPatch
 {
+    // Private data
+
+        //- Mapping engines
+        PtrList<directMappedPatchBase> mappers_;
+
 
     // Private Member Functions
 

src/sampling/sampledSurface/sampledPatchInternalField/sampledPatchInternalFieldTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,20 +36,27 @@ Foam::sampledPatchInternalField::sampleField
     const GeometricField<Type, fvPatchField, volMesh>& vField
 ) const
 {
-    const mapDistribute& distMap = map();
-
     // One value per face
     tmp<Field<Type> > tvalues(new Field<Type>(patchFaceLabels().size()));
     Field<Type>& values = tvalues();
 
-    if (patchIndex() != -1)
+    forAll(patchStart(), i)
     {
+        // Get patchface wise data by sampling internal field
         Field<Type> interpVals = vField.internalField();
-        distMap.distribute(interpVals);
+        mappers_[i].map().distribute(interpVals);
 
-        forAll(patchFaceLabels(), elemI)
+        // Store at correct position in values
+        label end =
+        (
+            i < patchStart().size()-1
+          ? patchStart()[i+1]
+          : patchFaceLabels().size()
+        );
+
+        for (label triI = patchStart()[i]; triI < end; triI++)
         {
-            values[elemI] = interpVals[patchFaceLabels()[elemI]];
+            values[triI] = interpVals[patchFaceLabels()[triI]];
         }
     }
 
@@ -64,19 +71,24 @@ Foam::sampledPatchInternalField::interpolateField
     const interpolation<Type>& interpolator
 ) const
 {
-    // One value per vertex
+    label sz = 0;
+    forAll(patchIDs(), i)
+    {
+        sz += mesh().boundaryMesh()[patchIDs()[i]].size();
+    }
 
-    if (patchIndex() != -1)
+    Field<Type> allPatchVals(sz);
+    sz = 0;
+
+    forAll(patchIDs(), i)
     {
         // See directMappedFixedValueFvPatchField
-        const mapDistribute& distMap = map();
-
-        const polyPatch& pp = mesh().boundaryMesh()[patchIndex()];
+        const mapDistribute& distMap = mappers_[i].map();
 
         // Send back sample points to processor that holds the cell.
         // Mark cells with point::max so we know which ones we need
         // to interpolate (since expensive).
-        vectorField samples(samplePoints());
+        vectorField samples(mappers_[i].samplePoints());
         distMap.reverseDistribute(mesh().nCells(), point::max, samples);
 
         Field<Type> patchVals(mesh().nCells());
@@ -96,18 +108,35 @@ Foam::sampledPatchInternalField::interpolateField
         distMap.distribute(patchVals);
 
         // Now patchVals holds the interpolated data in patch face order.
-        // Interpolate to points. Note: points are patch.localPoints() so
-        // can use standard interpolation
+        // Collect.
+        SubList<Type>(allPatchVals, patchVals.size(), sz).assign(patchVals);
+        sz += patchVals.size();
+    }
 
-        return PrimitivePatchInterpolation<primitivePatch>
-        (
-           pp
-        ).faceToPointInterpolate(patchVals);
-    }
-    else
+    // Interpolate to points. Reconstruct the patch of all faces to aid
+    // interpolation.
+
+    labelList meshFaceLabels(allPatchVals.size());
+    sz = 0;
+    forAll(patchIDs(), i)
     {
-        return tmp<Field<Type> >(new Field<Type>(points().size()));
+        const polyPatch& pp = mesh().boundaryMesh()[patchIDs()[i]];
+        forAll(pp, i)
+        {
+            meshFaceLabels[sz++] = pp.start()+i;
+        }
     }
+
+    indirectPrimitivePatch allPatches
+    (
+        IndirectList<face>(mesh().faces(), meshFaceLabels),
+        mesh().points()
+    );
+
+    return PrimitivePatchInterpolation<indirectPrimitivePatch>
+    (
+        allPatches
+    ).faceToPointInterpolate(allPatchVals);
 }
 
 

src/thermophysicalModels/radiationModels/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C

@@ -197,7 +197,7 @@ Foam::scalar Foam::radiation::radiativeIntensityRay::correct()
             1.0/constant::mathematical::pi*omega_
            *(
                 k*blackBody_.bLambda(lambdaI)
-              + absorptionEmission_.ECont(lambdaI)
+              + absorptionEmission_.ECont(lambdaI)/4
             )
         );
 

src/thermophysicalModels/thermalPorousZone/thermalModel/fixedTemperature/fixedTemperature.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -51,8 +51,8 @@ namespace porousMedia
 
 Foam::porousMedia::fixedTemperature::fixedTemperature(const porousZone& pZone)
 :
-    thermalModel(pZone, typeName),
-    T_(readScalar(coeffDict_.lookup("T")))
+    thermalModel(pZone),
+    T_(readScalar(thermalCoeffs_.lookup("T")))
 {}
 
 

src/thermophysicalModels/thermalPorousZone/thermalModel/fixedTemperature/fixedTemperature.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,7 +33,6 @@ Description
 #define fixedTemperature_H
 
 #include "thermalModel.H"
-#include "autoPtr.H"
 #include "runTimeSelectionTables.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -81,7 +80,7 @@ public:
         //- Add the thermal source to the enthalpy equation
         virtual void addEnthalpySource
         (
-            const basicThermo& thermo,
+            const basicThermo&,
             const volScalarField& rho,
             fvScalarMatrix& hEqn
         ) const;

src/thermophysicalModels/thermalPorousZone/thermalModel/noThermalModel/noThermalModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,7 +33,6 @@ Description
 #define noThermalModel_H
 
 #include "thermalModel.H"
-#include "autoPtr.H"
 #include "runTimeSelectionTables.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -73,7 +72,7 @@ public:
         //- Add the thermal source to the enthalpy equation
         virtual void addEnthalpySource
         (
-            const basicThermo& thermo,
+            const basicThermo&,
             const volScalarField& rho,
             fvScalarMatrix& hEqn
         ) const;

src/thermophysicalModels/thermalPorousZone/thermalModel/thermalModel/thermalModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -42,18 +42,18 @@ namespace porousMedia
 Foam::porousMedia::thermalModel::thermalModel(const porousZone& pZone)
 :
     pZone_(pZone),
-    coeffDict_(dictionary::null)
+    thermalCoeffs_(pZone.dict().subDictPtr("thermalModel"))
 {}
 
 
 Foam::porousMedia::thermalModel::thermalModel
 (
     const porousZone& pZone,
-    const word& modelType
+    const dictionary& thermalCoeffs
 )
 :
     pZone_(pZone),
-    coeffDict_(pZone_.dict().subDict(modelType + "Coeffs"))
+    thermalCoeffs_(thermalCoeffs)
 {}
 
 

src/thermophysicalModels/thermalPorousZone/thermalModel/thermalModel/thermalModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -63,8 +63,8 @@ protected:
         //- Reference to the porous zone
         const porousZone& pZone_;
 
-        //- Sub-model coefficients dictionary
-        const dictionary coeffDict_;
+        //- Thermal model coefficients dictionary
+        const dictionary thermalCoeffs_;
 
 
 public:
@@ -87,11 +87,11 @@ public:
 
     // Constructors
 
-        //- Construct null from porous zone
+        //- Construct from porous zone, coefficients from "thermalModel" entry
         thermalModel(const porousZone&);
 
-        //- Construct from porous zone and model type name
-        thermalModel(const porousZone&, const word& modelType);
+        //- Construct from porous zone and thermal model coefficients
+        thermalModel(const porousZone&, const dictionary& thermalCoeffs);
 
 
     //- Destructor
@@ -107,7 +107,7 @@ public:
         //- Add the thermal source to the enthalpy equation
         virtual void addEnthalpySource
         (
-            const basicThermo& thermo,
+            const basicThermo&,
             const volScalarField& rho,
             fvScalarMatrix& hEqn
         ) const = 0;

src/thermophysicalModels/thermalPorousZone/thermalModel/thermalModel/thermalModelNew.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,7 +33,13 @@ Foam::porousMedia::thermalModel::New
     const porousZone& pZone
 )
 {
-    const word modelType(pZone.dict().lookup("thermalModel"));
+    // a missing thermalModel is the same as type "none"
+    word modelType("none");
+
+    if (const dictionary* dictPtr = pZone.dict().subDictPtr("thermalModel"))
+    {
+        dictPtr->lookup("type") >> modelType;
+    }
 
     Info<< "Selecting thermalModel " << modelType << endl;
 

src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -51,7 +51,10 @@ void Foam::thermalPorousZone::addEnthalpySource
     fvScalarMatrix& hEqn
 ) const
 {
-    model_->addEnthalpySource(thermo, rho, hEqn);
+    if (model_.valid())
+    {
+        model_->addEnthalpySource(thermo, rho, hEqn);
+    }
 }
 
 

src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -127,7 +127,7 @@ public:
         //- Add the thermal source to the enthalpy equation
         void addEnthalpySource
         (
-            const basicThermo& thermo,
+            const basicThermo&,
             const volScalarField& rho,
             fvScalarMatrix& hEqn
         ) const;

src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZones.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZones.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -32,7 +32,8 @@ Description
     thermalPorousZone.
 
     The input file \c constant/thermalPorousZone is implemented as
-    IOPtrList\<thermalPorousZone\> and contains the following type of data:
+    IOPtrList\<thermalPorousZone\> (but written as a dictionary)
+    and contains the following type of data:
 
     \verbatim
     1
@@ -94,7 +95,7 @@ public:
         //- Add the thermal source to the enthalpy equation
         void addEnthalpySource
         (
-            const basicThermo& thermo,
+            const basicThermo&,
             const volScalarField& rho,
             fvScalarMatrix& hEqn
         ) const;

tutorials/basic/potentialFoam/cylinder/system/fvSchemes

@@ -22,7 +22,7 @@ ddtSchemes
 
 gradSchemes
 {
-    default         Gauss linear;
+    default         leastSquares;
 }
 
 divSchemes

tutorials/basic/potentialFoam/cylinder/system/fvSolution

@@ -26,7 +26,7 @@ solvers
     }
 }
 
-SIMPLE
+potentialFlow
 {
     nNonOrthogonalCorrectors 3;
 }

tutorials/basic/potentialFoam/pitzDaily/system/fvSolution

@@ -26,7 +26,7 @@ solvers
     }
 }
 
-SIMPLE
+potentialFlow
 {
     nNonOrthogonalCorrectors 0;
 }

tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/porousZones

@@ -31,11 +31,12 @@ FoamFile
             f   f [0 -1 0 0 0 0 0] (0 0 0);
         }
 
-        thermalModel none; // fixedTemperature;
-
-        fixedTemperatureCoeffs
+        thermalModel
         {
-            T              350;
+            type        none; // fixedTemperature;
+
+            // fixedTemperature coefficients
+            T           350;
         }
     }
 )

tutorials/incompressible/pisoFoam/les/pitzDaily/system/controlDict

@@ -124,7 +124,7 @@ functions
             nearWall
             {
                 type            patchInternalField;
-                patchName       lowerWall;
+                patches         ( lowerWall );
                 distance        1E-6;
                 interpolate     true;
                 triangulate     false;

tutorials/incompressible/pisoFoam/les/pitzDailyDirectMapped/system/controlDict

@@ -124,7 +124,7 @@ functions
             nearWall
             {
                 type            patchInternalField;
-                patchName       lowerWall;
+                patches         ( lowerWall );
                 distance        1E-6;
                 interpolate     true;
                 triangulate     false;

tutorials/incompressible/simpleFoam/motorBike/Allclean

@@ -3,7 +3,6 @@
 # Source tutorial clean functions
 . $WM_PROJECT_DIR/bin/tools/CleanFunctions
 
-rm -f system/fvSolution
 rm -rf 0 > /dev/null 2>&1
 
 cleanCase

tutorials/incompressible/simpleFoam/motorBike/Allrun

@@ -2,14 +2,10 @@
 # Source tutorial run functions
 . $WM_PROJECT_DIR/bin/tools/RunFunctions
 
-cp system/fvSolution.org system/fvSolution
 cp -r 0.org 0 > /dev/null 2>&1
 
 runApplication blockMesh
 runApplication snappyHexMesh -overwrite
 
-sed -i 's/\(nNonOrthogonalCorrectors\).*;/\1 10;/g' system/fvSolution
 runApplication potentialFoam -writep
-sed -i 's/\(nNonOrthogonalCorrectors\).*;/\1 0;/g' system/fvSolution
-
 runApplication `getApplication`

tutorials/incompressible/simpleFoam/motorBike/constant/polyMesh/boundary

@@ -21,433 +21,433 @@ FoamFile
     {
         type            patch;
         nFaces          320;
-        startFace       1028304;
+        startFace       1027872;
     }
     inlet
     {
         type            patch;
         nFaces          64;
-        startFace       1028624;
+        startFace       1028192;
     }
     outlet
     {
         type            patch;
         nFaces          64;
-        startFace       1028688;
+        startFace       1028256;
     }
     lowerWall
     {
         type            wall;
         nFaces          5330;
-        startFace       1028752;
+        startFace       1028320;
     }
     upperWall
     {
         type            patch;
         nFaces          160;
-        startFace       1034082;
+        startFace       1033650;
     }
     motorBike_frt-fairing:001%1
     {
         type            wall;
-        nFaces          5427;
-        startFace       1034242;
+        nFaces          5416;
+        startFace       1033810;
     }
     motorBike_windshield:002%2
     {
         type            wall;
         nFaces          50;
-        startFace       1039669;
+        startFace       1039226;
     }
     motorBike_rr-wh-rim:005%5
     {
         type            wall;
-        nFaces          136;
-        startFace       1039719;
+        nFaces          137;
+        startFace       1039276;
     }
     motorBike_rr-wh-rim:010%10
     {
         type            wall;
         nFaces          340;
-        startFace       1039855;
+        startFace       1039413;
     }
     motorBike_fr-wh-rim:011%11
     {
         type            wall;
         nFaces          473;
-        startFace       1040195;
+        startFace       1039753;
     }
     motorBike_fr-wh-brake-disk:012%12
     {
         type            wall;
         nFaces          54;
-        startFace       1040668;
+        startFace       1040226;
     }
     motorBike_frame:016-shadow%13
     {
         type            wall;
-        nFaces          97;
-        startFace       1040722;
+        nFaces          98;
+        startFace       1040280;
     }
     motorBike_rear-susp:014%14
     {
         type            wall;
-        nFaces          839;
-        startFace       1040819;
+        nFaces          844;
+        startFace       1040378;
     }
     motorBike_rear-susp:014-shadow%15
     {
         type            wall;
-        nFaces          469;
-        startFace       1041658;
+        nFaces          461;
+        startFace       1041222;
     }
     motorBike_frame:016%16
     {
         type            wall;
-        nFaces          68;
-        startFace       1042127;
+        nFaces          67;
+        startFace       1041683;
     }
     motorBike_rr-wh-rim:005-shadow%17
     {
         type            wall;
-        nFaces          71;
-        startFace       1042195;
+        nFaces          70;
+        startFace       1041750;
     }
     motorBike_rr-wh-chain-hub:022%22
     {
         type            wall;
         nFaces          141;
-        startFace       1042266;
+        startFace       1041820;
     }
     motorBike_rearseat%24
     {
         type            wall;
         nFaces          430;
-        startFace       1042407;
+        startFace       1041961;
     }
     motorBike_frt-fairing%25
     {
         type            wall;
         nFaces          626;
-        startFace       1042837;
+        startFace       1042391;
     }
     motorBike_windshield%26
     {
         type            wall;
-        nFaces          378;
-        startFace       1043463;
+        nFaces          368;
+        startFace       1043017;
     }
     motorBike_headlights%27
     {
         type            wall;
         nFaces          161;
-        startFace       1043841;
+        startFace       1043385;
     }
     motorBike_driversseat%28
     {
         type            wall;
         nFaces          368;
-        startFace       1044002;
+        startFace       1043546;
     }
     motorBike_rear-body%29
     {
         type            wall;
         nFaces          2077;
-        startFace       1044370;
+        startFace       1043914;
     }
     motorBike_fuel-tank%30
     {
         type            wall;
         nFaces          912;
-        startFace       1046447;
+        startFace       1045991;
     }
     motorBike_exhaust%31
     {
         type            wall;
         nFaces          2392;
-        startFace       1047359;
+        startFace       1046903;
     }
     motorBike_rr-wh-rim%32
     {
         type            wall;
         nFaces          1430;
-        startFace       1049751;
+        startFace       1049295;
     }
     motorBike_fr-mud-guard%33
     {
         type            wall;
-        nFaces          639;
-        startFace       1051181;
+        nFaces          624;
+        startFace       1050725;
     }
     motorBike_fr-wh-rim%34
     {
         type            wall;
         nFaces          591;
-        startFace       1051820;
+        startFace       1051349;
     }
     motorBike_fr-wh-brake-disk%35
     {
         type            wall;
-        nFaces          431;
-        startFace       1052411;
+        nFaces          427;
+        startFace       1051940;
     }
     motorBike_fr-brake-caliper%36
     {
         type            wall;
         nFaces          164;
-        startFace       1052842;
+        startFace       1052367;
     }
     motorBike_fr-wh-tyre%37
     {
         type            wall;
         nFaces          1116;
-        startFace       1053006;
+        startFace       1052531;
     }
     motorBike_hbars%38
     {
         type            wall;
         nFaces          535;
-        startFace       1054122;
+        startFace       1053647;
     }
     motorBike_fr-forks%39
     {
         type            wall;
         nFaces          1140;
-        startFace       1054657;
+        startFace       1054182;
     }
     motorBike_chain%40
     {
         type            wall;
         nFaces          474;
-        startFace       1055797;
+        startFace       1055322;
     }
     motorBike_rr-wh-tyre%41
     {
         type            wall;
         nFaces          1787;
-        startFace       1056271;
+        startFace       1055796;
     }
     motorBike_square-dial%42
     {
         type            wall;
         nFaces          6;
-        startFace       1058058;
+        startFace       1057583;
     }
     motorBike_round-dial%43
     {
         type            wall;
         nFaces          17;
-        startFace       1058064;
+        startFace       1057589;
     }
     motorBike_dial-holder%44
     {
         type            wall;
         nFaces          87;
-        startFace       1058081;
+        startFace       1057606;
     }
     motorBike_rear-susp%45
     {
         type            wall;
         nFaces          1786;
-        startFace       1058168;
+        startFace       1057693;
     }
     motorBike_rear-brake-lights%46
     {
         type            wall;
         nFaces          53;
-        startFace       1059954;
+        startFace       1059479;
     }
     motorBike_rear-light-bracket%47
     {
         type            wall;
         nFaces          163;
-        startFace       1060007;
+        startFace       1059532;
     }
     motorBike_frame%48
     {
         type            wall;
         nFaces          2040;
-        startFace       1060170;
+        startFace       1059695;
     }
     motorBike_rear-mud-guard%49
     {
         type            wall;
         nFaces          663;
-        startFace       1062210;
+        startFace       1061735;
     }
     motorBike_rear-susp-spring-damp%50
     {
         type            wall;
         nFaces          107;
-        startFace       1062873;
+        startFace       1062398;
     }
     motorBike_fairing-inner-plate%51
     {
         type            wall;
         nFaces          445;
-        startFace       1062980;
+        startFace       1062505;
     }
     motorBike_clutch-housing%52
     {
         type            wall;
         nFaces          970;
-        startFace       1063425;
+        startFace       1062950;
     }
     motorBike_radiator%53
     {
         type            wall;
-        nFaces          40;
-        startFace       1064395;
+        nFaces          39;
+        startFace       1063920;
     }
     motorBike_water-pipe%54
     {
         type            wall;
         nFaces          104;
-        startFace       1064435;
+        startFace       1063959;
     }
     motorBike_water-pump%55
     {
         type            wall;
         nFaces          74;
-        startFace       1064539;
+        startFace       1064063;
     }
     motorBike_engine%56
     {
         type            wall;
         nFaces          2389;
-        startFace       1064613;
+        startFace       1064137;
     }
     motorBike_rear-shock-link%57
     {
         type            wall;
-        nFaces          25;
-        startFace       1067002;
+        nFaces          26;
+        startFace       1066526;
     }
     motorBike_rear-brake-fluid-pot-bracket%58
     {
         type            wall;
         nFaces          44;
-        startFace       1067027;
+        startFace       1066552;
     }
     motorBike_rear-brake-fluid-pot%59
     {
         type            wall;
         nFaces          53;
-        startFace       1067071;
+        startFace       1066596;
     }
     motorBike_footpeg%60
     {
         type            wall;
         nFaces          86;
-        startFace       1067124;
+        startFace       1066649;
     }
     motorBike_rr-wh-chain-hub%61
     {
         type            wall;
         nFaces          122;
-        startFace       1067210;
+        startFace       1066735;
     }
     motorBike_rear-brake-caliper%62
     {
         type            wall;
         nFaces          142;
-        startFace       1067332;
+        startFace       1066857;
     }
     motorBike_rider-helmet%65
     {
         type            wall;
         nFaces          583;
-        startFace       1067474;
+        startFace       1066999;
     }
     motorBike_rider-visor%66
     {
         type            wall;
         nFaces          95;
-        startFace       1068057;
+        startFace       1067582;
     }
     motorBike_rider-boots%67
     {
         type            wall;
         nFaces          1025;
-        startFace       1068152;
+        startFace       1067677;
     }
     motorBike_rider-gloves%68
     {
         type            wall;
         nFaces          319;
-        startFace       1069177;
+        startFace       1068702;
     }
     motorBike_rider-body%69
     {
         type            wall;
         nFaces          4555;
-        startFace       1069496;
+        startFace       1069021;
     }
     motorBike_frame:0%70
     {
         type            wall;
         nFaces          37;
-        startFace       1074051;
+        startFace       1073576;
     }
     motorBike_frt-fairing:001-shadow%74
     {
         type            wall;
-        nFaces          3306;
-        startFace       1074088;
+        nFaces          3317;
+        startFace       1073613;
     }
     motorBike_windshield-shadow%75
     {
         type            wall;
-        nFaces          239;
-        startFace       1077394;
+        nFaces          236;
+        startFace       1076930;
     }
     motorBike_fr-mud-guard-shadow%81
     {
         type            wall;
-        nFaces          332;
-        startFace       1077633;
+        nFaces          350;
+        startFace       1077166;
     }
     motorBike_fr-wh-brake-disk-shadow%83
     {
         type            wall;
-        nFaces          182;
-        startFace       1077965;
+        nFaces          186;
+        startFace       1077516;
     }
     motorBike_rear-mud-guard-shadow%84
     {
         type            wall;
-        nFaces          387;
-        startFace       1078147;
+        nFaces          388;
+        startFace       1077702;
     }
     motorBike_rear-susp-spring-damp-shadow%85
     {
         type            wall;
         nFaces          41;
-        startFace       1078534;
+        startFace       1078090;
     }
     motorBike_radiator-shadow%86
     {
         type            wall;
-        nFaces          21;
-        startFace       1078575;
+        nFaces          22;
+        startFace       1078131;
     }
     motorBike_rear-shock-link-shadow%87
     {
         type            wall;
-        nFaces          10;
-        startFace       1078596;
+        nFaces          9;
+        startFace       1078153;
     }
     motorBike_rear-brake-fluid-pot-bracket-shadow%88
     {
         type            wall;
         nFaces          22;
-        startFace       1078606;
+        startFace       1078162;
     }
     motorBike_rr-wh-chain-hub-shadow%89
     {
         type            wall;
         nFaces          53;
-        startFace       1078628;
+        startFace       1078184;
     }
 )
 

tutorials/incompressible/simpleFoam/motorBike/system/fvSolution

@@ -63,6 +63,11 @@ SIMPLE
     nNonOrthogonalCorrectors 0;
 }
 
+potentialFlow
+{
+    nNonOrthogonalCorrectors 10;
+}
+
 relaxationFactors
 {
     p               0.3;

tutorials/incompressible/simpleFoam/motorBike/system/snappyHexMeshDict

@@ -175,7 +175,7 @@ snapControls
     tolerance 4.0;
 
     //- Number of mesh displacement relaxation iterations.
-    nSolveIter 30;
+    nSolveIter 0;
 
     //- Maximum number of snapping relaxation iterations. Should stop
     //  before upon reaching a correct mesh.

tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSolution

@@ -74,7 +74,7 @@ PISO
     momentumPredictor yes;
 }
 
-SIMPLE
+potentialFlow
 {
     // used for potentialFoam initialisation
     nNonOrthogonalCorrectors 10;

tutorials/mesh/snappyHexMesh/flange/Allclean

@@ -0,0 +1,15 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+rm -rf 0 > /dev/null 2>&1
+rm -f ./*.obj > /dev/null 2>&1
+rm -rf constant/extendedFeatureEdgeMesh > /dev/null 2>&1
+rm -f constant/triSurface/flange.eMesh > /dev/null 2>&1
+rm -f constant/polyMesh/boundary
+
+cleanCase
+
+# ----------------------------------------------------------------- end-of-file

tutorials/mesh/snappyHexMesh/flange/Allrun

@@ -0,0 +1,7 @@
+#!/bin/sh
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+runApplication blockMesh
+runApplication surfaceFeatureExtract -includedAngle 135 -writeObj constant/triSurface/flange.stl flange
+runApplication snappyHexMesh -overwrite

tutorials/mesh/snappyHexMesh/flange/constant/polyMesh/blockMeshDict

@@ -0,0 +1,53 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 1;
+
+vertices
+(
+    ( -0.04 -0.043 -0.04)
+    (  0.04 -0.043 -0.04)
+    (  0.04  0.04 -0.04)
+    ( -0.04  0.04 -0.04)
+    ( -0.04 -0.043 0.01)
+    (  0.04 -0.043 0.01)
+    (  0.04  0.04 0.01)
+    ( -0.04  0.04 0.01)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (20 20 12) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+patches
+(
+    patch allBoundary
+    (
+        (3 7 6 2)
+        (0 4 7 3)
+        (2 6 5 1)
+        (1 5 4 0)
+        (0 3 2 1)
+        (4 5 6 7)
+    )
+);
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/flange/constant/transportProperties

@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+nu              nu [ 0 2 -1 0 0 0 0 ] 0.01;
+
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/flange/system/controlDict

@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     icoFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         50;
+
+deltaT          1;
+
+writeControl    timeStep;
+
+writeInterval   20;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/flange/system/decomposeParDict

@@ -0,0 +1,50 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 4;
+
+method          ptscotch; //simple;
+
+simpleCoeffs
+{
+    n               ( 2 2 1 );
+    delta           0.001;
+}
+
+hierarchicalCoeffs
+{
+    n               ( 1 1 1 );
+    delta           0.001;
+    order           xyz;
+}
+
+metisCoeffs
+{
+    processorWeights ( 1 1 1 1 );
+}
+
+manualCoeffs
+{
+    dataFile        "";
+}
+
+distributed     no;
+
+roots           ( );
+
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/flange/system/fvSchemes

@@ -0,0 +1,60 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+    grad(p)         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+    div(phi,U)      Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         none;
+    laplacian(nu,U) Gauss linear corrected;
+    laplacian((1|A(U)),p) Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+    interpolate(HbyA) linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p               ;
+}
+
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/flange/system/fvSolution

@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-06;
+        relTol          0;
+    }
+
+    U
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-05;
+        relTol          0;
+    }
+}
+
+PISO
+{
+    nCorrectors     2;
+    nNonOrthogonalCorrectors 0;
+    pRefCell        0;
+    pRefValue       0;
+}
+
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/flange/system/snappyHexMeshDict

@@ -0,0 +1,357 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                  |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      autoHexMeshDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Which of the steps to run
+castellatedMesh true;
+snap            true;
+addLayers       false;
+
+
+// Geometry. Definition of all surfaces. All surfaces are of class
+// searchableSurface.
+// Surfaces are used
+// - to specify refinement for any mesh cell intersecting it
+// - to specify refinement for any mesh cell inside/outside/near
+// - to 'snap' the mesh boundary to the surface
+geometry
+{
+    flange.stl
+    {
+        type triSurfaceMesh;
+        name flange;
+    }
+
+    refine
+    {
+        type searchableSphere;
+        centre (0 0 -0.012);
+        radius 0.003;
+    }
+};
+
+
+
+// Settings for the castellatedMesh generation.
+castellatedMeshControls
+{
+
+    // Refinement parameters
+    // ~~~~~~~~~~~~~~~~~~~~~
+
+    // While refining maximum number of cells per processor. This is basically
+    // the number of cells that fit on a processor. If you choose this too small
+    // it will do just more refinement iterations to obtain a similar mesh.
+    maxLocalCells 1000000;
+
+    // Overall cell limit (approximately). Refinement will stop immediately
+    // upon reaching this number so a refinement level might not complete.
+    // Note that this is the number of cells before removing the part which
+    // is not 'visible' from the keepPoint. The final number of cells might
+    // actually be a lot less.
+    maxGlobalCells 2000000;
+
+    // The surface refinement loop might spend lots of iterations refining just a
+    // few cells. This setting will cause refinement to stop if <= minimumRefine
+    // are selected for refinement. Note: it will at least do one iteration
+    // (unless the number of cells to refine is 0)
+    minRefinementCells 0;
+
+    // Number of buffer layers between different levels.
+    // 1 means normal 2:1 refinement restriction, larger means slower
+    // refinement.
+    nCellsBetweenLevels 1;
+
+
+
+    // Explicit feature edge refinement
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Specifies a level for any cell intersected by its edges.
+    // This is a featureEdgeMesh, read from constant/triSurface for now.
+    features
+    (
+        {
+            file "flange.eMesh";
+            level 0;
+        }
+    );
+
+
+
+    // Surface based refinement
+    // ~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Specifies two levels for every surface. The first is the minimum level,
+    // every cell intersecting a surface gets refined up to the minimum level.
+    // The second level is the maximum level. Cells that 'see' multiple
+    // intersections where the intersections make an
+    // angle > resolveFeatureAngle get refined up to the maximum level.
+
+    refinementSurfaces
+    {
+        flange
+        {
+            // Surface-wise min and max refinement level
+            level (2 2);
+        }
+    }
+
+    resolveFeatureAngle 45;
+
+
+    // Region-wise refinement
+    // ~~~~~~~~~~~~~~~~~~~~~~
+
+    // Specifies refinement level for cells in relation to a surface. One of
+    // three modes
+    // - distance. 'levels' specifies per distance to the surface the
+    //   wanted refinement level. The distances need to be specified in
+    //   descending order.
+    // - inside. 'levels' is only one entry and only the level is used. All
+    //   cells inside the surface get refined up to the level. The surface
+    //   needs to be closed for this to be possible.
+    // - outside. Same but cells outside.
+
+    refinementRegions
+    {
+        refine
+        {
+            mode inside;
+            levels ((1E15 3));
+        }
+    }
+
+
+    // Mesh selection
+    // ~~~~~~~~~~~~~~
+
+    // 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.
+    // NOTE: This point should never be on a face, always inside a cell, even
+    // after refinement.
+    // This is an outside point locationInMesh (-0.033 -0.033 0.0033);
+    locationInMesh (-9.23149e-05 -0.0025 -0.0025); // Inside point
+
+    // 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;
+}
+
+
+
+// Settings for the snapping.
+snapControls
+{
+    //- Number of patch smoothing iterations before finding correspondence
+    //  to surface
+    nSmoothPatch 3;
+
+    //- Relative distance for points to be attracted by surface feature point
+    //  or edge. True distance is this factor times local
+    //  maximum edge length.
+    tolerance 1.0;
+
+    //- Number of mesh displacement relaxation iterations.
+    nSolveIter 300;
+
+    //- Maximum number of snapping relaxation iterations. Should stop
+    //  before upon reaching a correct mesh.
+    nRelaxIter 5;
+
+    //- Highly experimental and wip: number of feature edge snapping
+    //  iterations. Leave out altogether to disable.
+    nFeatureSnapIter 10;
+}
+
+
+
+// Settings for the layer addition.
+addLayersControls
+{
+    // Are the thickness parameters below relative to the undistorted
+    // size of the refined cell outside layer (true) or absolute sizes (false).
+    relativeSizes true;
+
+    // Per final patch (so not geometry!) the layer information
+    layers
+    {
+        "flange_.*"
+        {
+            nSurfaceLayers 1;
+        }
+    }
+
+    // Expansion factor for layer mesh
+    expansionRatio 1.0;
+
+
+    //- Wanted thickness of final added cell layer. If multiple layers
+    //  is the thickness of the layer furthest away from the wall.
+    //  See relativeSizes parameter.
+    finalLayerThickness 0.3;
+
+    //- Minimum thickness of cell layer. If for any reason layer
+    //  cannot be above minThickness do not add layer.
+    //  See relativeSizes parameter.
+    minThickness 0.25;
+
+    //- If points get not extruded do nGrow layers of connected faces that are
+    //  also not grown. This helps convergence of the layer addition process
+    //  close to features.
+    nGrow 0;
+
+
+    // Advanced settings
+
+    //- When not to extrude surface. 0 is flat surface, 90 is when two faces
+    //  make straight angle.
+    featureAngle 45;
+
+    //- Maximum number of snapping relaxation iterations. Should stop
+    //  before upon reaching a correct mesh.
+    nRelaxIter 5;
+
+    // Number of smoothing iterations of surface normals
+    nSmoothSurfaceNormals 1;
+
+    // Number of smoothing iterations of interior mesh movement direction
+    nSmoothNormals 3;
+
+    // Smooth layer thickness over surface patches
+    nSmoothThickness 10;
+
+    // Stop layer growth on highly warped cells
+    maxFaceThicknessRatio 0.5;
+
+    // Reduce layer growth where ratio thickness to medial
+    // distance is large
+    maxThicknessToMedialRatio 0.3;
+
+    // Angle used to pick up medial axis points
+    minMedianAxisAngle 90;
+
+    // Create buffer region for new layer terminations
+    nBufferCellsNoExtrude 0;
+
+
+    // Overall max number of layer addition iterations. The mesher will exit
+    // if it reaches this number of iterations; possibly with an illegal
+    // mesh.
+    nLayerIter 50;
+
+    // Max number of iterations after which relaxed meshQuality controls
+    // get used. Up to nRelaxIter it uses the settings in meshQualityControls,
+    // after nRelaxIter it uses the values in meshQualityControls::relaxed.
+    nRelaxedIter 20;
+}
+
+
+
+// Generic mesh quality settings. At any undoable phase these determine
+// where to undo.
+meshQualityControls
+{
+    //- Maximum non-orthogonality allowed. Set to 180 to disable.
+    maxNonOrtho 65;
+
+    //- Max skewness allowed. Set to <0 to disable.
+    maxBoundarySkewness 20;
+    maxInternalSkewness 4;
+
+    //- Max concaveness allowed. Is angle (in degrees) below which concavity
+    //  is allowed. 0 is straight face, <0 would be convex face.
+    //  Set to 180 to disable.
+    maxConcave 80;
+
+    //- Minimum pyramid volume. Is absolute volume of cell pyramid.
+    //  Set to a sensible fraction of the smallest cell volume expected.
+    //  Set to very negative number (e.g. -1E30) to disable.
+    minVol 1e-13;
+
+    //- Minimum quality of the tet formed by the face-centre
+    //  and variable base point minimum decomposition triangles and
+    //  the cell centre.  Set to very negative number (e.g. -1E30) to
+    //  disable.
+    //     <0 = inside out tet,
+    //      0 = flat tet
+    //      1 = regular tet
+    minTetQuality 1e-9;
+
+    //- Minimum face area. Set to <0 to disable.
+    minArea -1;
+
+    //- Minimum face twist. Set to <-1 to disable. dot product of face normal
+    //- and face centre triangles normal
+    minTwist 0.05;
+
+    //- minimum normalised cell determinant
+    //- 1 = hex, <= 0 = folded or flattened illegal cell
+    minDeterminant 0.001;
+
+    //- minFaceWeight (0 -> 0.5)
+    minFaceWeight 0.05;
+
+    //- minVolRatio (0 -> 1)
+    minVolRatio 0.01;
+
+    //must be >0 for Fluent compatibility
+    minTriangleTwist -1;
+
+    //- if >0 : preserve single cells with all points on the surface if the
+    //  resulting volume after snapping (by approximation) is larger than
+    //  minVolCollapseRatio times old volume (i.e. not collapsed to flat cell).
+    //  If <0 : delete always.
+    //minVolCollapseRatio 0.5;
+
+
+    // Advanced
+
+    //- Number of error distribution iterations
+    nSmoothScale 4;
+    //- amount to scale back displacement at error points
+    errorReduction 0.75;
+
+
+
+    // Optional : some meshing phases allow usage of relaxed rules.
+    // See e.g. addLayersControls::nRelaxedIter.
+    relaxed
+    {
+        //- Maximum non-orthogonality allowed. Set to 180 to disable.
+        maxNonOrtho 75;
+    }
+}
+
+
+// Advanced
+
+// Flags for optional output
+// 0 : only write final meshes
+// 1 : write intermediate meshes
+// 2 : write volScalarField with cellLevel for postprocessing
+// 4 : write current intersections as .obj files
+debug 0;
+
+
+// Merge tolerance. Is fraction of overall bounding box of initial mesh.
+// Note: the write tolerance needs to be higher than this.
+mergeTolerance 1E-6;
+
+
+// ************************************************************************* //

tutorials/multiphase/interDyMFoam/ras/sloshingTank2D/system/controlDict

@@ -89,7 +89,7 @@ functions
             walls
             {
                 type        patch;
-                patchName   walls;
+                patches     (walls);
                 triangulate false;
             }
         );

tutorials/multiphase/interDyMFoam/ras/sloshingTank2D3DoF/system/controlDict

@@ -88,7 +88,7 @@ functions
             walls
             {
                 type        patch;
-                patchName   walls;
+                patches     (walls);
                 triangulate false;
             }
         );

wmake/rules/linux64Gcc46/c

@@ -0,0 +1,16 @@
+.SUFFIXES: .c .h
+
+cWARN        = -Wall
+
+cc          = gcc -m64
+
+include $(RULES)/c$(WM_COMPILE_OPTION)
+
+cFLAGS      = $(GFLAGS) $(cWARN) $(cOPT) $(cDBUG) $(LIB_HEADER_DIRS) -fPIC
+
+ctoo        = $(WM_SCHEDULER) $(cc) $(cFLAGS) -c $$SOURCE -o $@
+
+LINK_LIBS   = $(cDBUG)
+
+LINKLIBSO   = $(cc) -shared
+LINKEXE     = $(cc) -Xlinker --add-needed -Xlinker -z -Xlinker nodefs

wmake/rules/linux64Gcc46/c++

@@ -0,0 +1,21 @@
+.SUFFIXES: .C .cxx .cc .cpp
+
+c++WARN     = -Wall -Wextra -Wno-unused-parameter -Wold-style-cast
+
+CC          = g++ -m64
+
+include $(RULES)/c++$(WM_COMPILE_OPTION)
+
+ptFLAGS     = -DNoRepository -ftemplate-depth-100
+
+c++FLAGS    = $(GFLAGS) $(c++WARN) $(c++OPT) $(c++DBUG) $(ptFLAGS) $(LIB_HEADER_DIRS) -fPIC
+
+Ctoo        = $(WM_SCHEDULER) $(CC) $(c++FLAGS) -c $$SOURCE -o $@
+cxxtoo      = $(Ctoo)
+cctoo       = $(Ctoo)
+cpptoo      = $(Ctoo)
+
+LINK_LIBS   = $(c++DBUG)
+
+LINKLIBSO   = $(CC) $(c++FLAGS) -shared
+LINKEXE     = $(CC) $(c++FLAGS) -Xlinker --add-needed

wmake/rules/linux64Gcc46/c++Debug

@@ -0,0 +1,2 @@
+c++DBUG    = -ggdb3 -DFULLDEBUG
+c++OPT      = -O0 -fdefault-inline

wmake/rules/linux64Gcc46/c++Opt

@@ -0,0 +1,4 @@
+c++DBUG     =
+c++OPT      = -O3
+#c++OPT      = -march=nocona -O3
+# -ftree-vectorize -ftree-vectorizer-verbose=3

wmake/rules/linux64Gcc46/c++Prof

@@ -0,0 +1,2 @@
+c++DBUG    = -pg
+c++OPT     = -O2

wmake/rules/linux64Gcc46/cDebug

@@ -0,0 +1,2 @@
+cDBUG       = -ggdb -DFULLDEBUG
+cOPT        = -O1 -fdefault-inline -finline-functions

wmake/rules/linux64Gcc46/cOpt

@@ -0,0 +1,2 @@
+cDBUG       =
+cOPT        = -O3

wmake/rules/linux64Gcc46/cProf

@@ -0,0 +1,2 @@
+cDBUG       = -pg
+cOPT        = -O2

wmake/rules/linux64Gcc46/general

@@ -0,0 +1,8 @@
+CPP        = cpp -traditional-cpp
+
+PROJECT_LIBS = -l$(WM_PROJECT) -ldl
+
+include $(GENERAL_RULES)/standard
+
+include $(RULES)/c
+include $(RULES)/c++

wmake/rules/linux64Gcc46/mplibHPMPI

@@ -0,0 +1,3 @@
+PFLAGS     =
+PINC       = -I$(MPI_ARCH_PATH)/include -D_MPICC_H
+PLIBS      = -L$(MPI_ARCH_PATH)/lib/linux_amd64 -lmpi

wmake/rules/linuxGcc46/c

@@ -0,0 +1,16 @@
+.SUFFIXES: .c .h
+
+cWARN        = -Wall
+
+cc          = gcc -m32
+
+include $(RULES)/c$(WM_COMPILE_OPTION)
+
+cFLAGS      = $(GFLAGS) $(cWARN) $(cOPT) $(cDBUG) $(LIB_HEADER_DIRS) -fPIC
+
+ctoo        = $(WM_SCHEDULER) $(cc) $(cFLAGS) -c $$SOURCE -o $@
+
+LINK_LIBS   = $(cDBUG)
+
+LINKLIBSO   = $(cc) -shared
+LINKEXE     = $(cc) -Xlinker --add-needed -Xlinker -z -Xlinker nodefs

wmake/rules/linuxGcc46/c++

@@ -0,0 +1,21 @@
+.SUFFIXES: .C .cxx .cc .cpp
+
+c++WARN     = -Wall -Wextra -Wno-unused-parameter -Wold-style-cast
+
+CC          = g++ -m32
+
+include $(RULES)/c++$(WM_COMPILE_OPTION)
+
+ptFLAGS     = -DNoRepository -ftemplate-depth-100
+
+c++FLAGS    = $(GFLAGS) $(c++WARN) $(c++OPT) $(c++DBUG) $(ptFLAGS) $(LIB_HEADER_DIRS) -fPIC
+
+Ctoo        = $(WM_SCHEDULER) $(CC) $(c++FLAGS) -c $$SOURCE -o $@
+cxxtoo      = $(Ctoo)
+cctoo       = $(Ctoo)
+cpptoo      = $(Ctoo)
+
+LINK_LIBS   = $(c++DBUG)
+
+LINKLIBSO   = $(CC) $(c++FLAGS) -shared
+LINKEXE     = $(CC) $(c++FLAGS) -Xlinker --add-needed

wmake/rules/linuxGcc46/c++Debug

@@ -0,0 +1,2 @@
+c++DBUG    = -ggdb3 -DFULLDEBUG
+c++OPT      = -O0 -fdefault-inline

wmake/rules/linuxGcc46/c++Opt

@@ -0,0 +1,2 @@
+c++DBUG     =
+c++OPT      = -O3

wmake/rules/linuxGcc46/c++Prof

@@ -0,0 +1,2 @@
+c++DBUG    = -pg
+c++OPT     = -O2