diff --git a/applications/solvers/combustion/chemFoam/readInitialConditions.H b/applications/solvers/combustion/chemFoam/readInitialConditions.H
index 3fbb13a0fc..08050be816 100644
--- a/applications/solvers/combustion/chemFoam/readInitialConditions.H
+++ b/applications/solvers/combustion/chemFoam/readInitialConditions.H
@@ -99,7 +99,7 @@
     scalar u0 = hs0 - p0/rho0;
     scalar R0 = p0/(rho0*T0);
     Rspecific[0] = R0;
- 
+
     scalar integratedHeat = 0.0;
 
     Info << constProp << " will be held constant." << nl
diff --git a/applications/test/vectorTools/Test-vectorTools.C b/applications/test/vectorTools/Test-vectorTools.C
index 0c85ac935b..85d18ed989 100644
--- a/applications/test/vectorTools/Test-vectorTools.C
+++ b/applications/test/vectorTools/Test-vectorTools.C
@@ -8,7 +8,7 @@ using namespace Foam;
 
 void test(const vector& a, const vector& b, const scalar tolerance)
 {
-    Info<< "Vectors " << a << " and " << b 
+    Info<< "Vectors " << a << " and " << b
         << " are (to tolerance of " << tolerance << "): ";
 
     if (vectorTools::areParallel(a, b, tolerance))
diff --git a/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/cv2DControls/cv2DControls.H b/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/cv2DControls/cv2DControls.H
index 4ab8bdb167..0b3cf0b4f2 100644
--- a/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/cv2DControls/cv2DControls.H
+++ b/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/cv2DControls/cv2DControls.H
@@ -252,4 +252,4 @@ public:
 
 #endif
 
-// ************************************************************************* // 
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cv2DMesh/insertSurfaceNearPointPairs.C b/applications/utilities/mesh/generation/cv2DMesh/insertSurfaceNearPointPairs.C
index a1d1919f12..3d4ac7cf2a 100644
--- a/applications/utilities/mesh/generation/cv2DMesh/insertSurfaceNearPointPairs.C
+++ b/applications/utilities/mesh/generation/cv2DMesh/insertSurfaceNearPointPairs.C
@@ -93,7 +93,7 @@ void Foam::CV2D::insertSurfaceNearPointPairs()
 
                         nNearPoints++;
 
-                        // Correct the edge iterator for the change in the 
+                        // Correct the edge iterator for the change in the
                         // number of edges following the point-pair insertion
                         eit = Finite_edges_iterator
                         (
diff --git a/applications/utilities/mesh/generation/cvMesh/Make/options b/applications/utilities/mesh/generation/cvMesh/Make/options
index a0a776d977..5efdb12572 100644
--- a/applications/utilities/mesh/generation/cvMesh/Make/options
+++ b/applications/utilities/mesh/generation/cvMesh/Make/options
@@ -20,7 +20,8 @@ EXE_LIBS = \
     $(CGAL_LIBS) \
     -lconformalVoronoiMesh \
     -lmeshTools \
-    -ldecompositionMethods -L$(FOAM_LIBBIN)/dummy -lscotchDecomp \
+    -ldecompositionMethods \
+    -L$(FOAM_LIBBIN)/dummy -lptscotchDecomp \
     -ledgeMesh \
     -ltriSurface \
     -ldynamicMesh
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/options b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/options
index 1a1902b559..2d2d1a1178 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/options
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/options
@@ -19,7 +19,6 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lmeshTools \
-    -ldecompositionMethods -L$(FOAM_LIBBIN)/dummy -lscotchDecomp \
     -ledgeMesh \
     -ltriSurface \
     -ldynamicMesh
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/cellSizeCalculationType/fieldFromFile/fieldFromFile.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/cellSizeCalculationType/fieldFromFile/fieldFromFile.C
index d40427abf7..788c1d69c6 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/cellSizeCalculationType/fieldFromFile/fieldFromFile.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/cellSizeCalculationType/fieldFromFile/fieldFromFile.C
@@ -63,7 +63,7 @@ Foam::fieldFromFile::fieldFromFile
 
 Foam::triSurfaceScalarField Foam::fieldFromFile::load()
 {
-    Info<< "Loading: " << fileName_ << endl;
+    Info<< indent << "Loading: " << fileName_ << endl;
 
     triSurfaceScalarField surfaceCellSize
     (
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/nonUniformField/nonUniformField.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/nonUniformField/nonUniformField.C
index 5e2d032d91..fc91ed279b 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/nonUniformField/nonUniformField.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/nonUniformField/nonUniformField.C
@@ -82,9 +82,9 @@ Foam::nonUniformField::nonUniformField
     Info<< decrIndent;
 
     Info<< indent << "Cell size field statistics:" << nl
-        << indent << "   Minimum: " << min(surfaceCellSize_).value() << nl
-        << indent << "   Average: " << average(surfaceCellSize_).value() << nl
-        << indent << "   Maximum: " << max(surfaceCellSize_).value() << endl;
+        << indent << "    Minimum: " << min(surfaceCellSize_).value() << nl
+        << indent << "    Average: " << average(surfaceCellSize_).value() << nl
+        << indent << "    Maximum: " << max(surfaceCellSize_).value() << endl;
 
     Info<< decrIndent;
 }
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/surfaceCellSizeFunction/surfaceCellSizeFunction.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/surfaceCellSizeFunction/surfaceCellSizeFunction.C
index 7022511c4c..04d9af824f 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/surfaceCellSizeFunction/surfaceCellSizeFunction.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/surfaceCellSizeFunction/surfaceCellSizeFunction.C
@@ -46,7 +46,10 @@ Foam::surfaceCellSizeFunction::surfaceCellSizeFunction
     dictionary(surfaceCellSizeFunctionDict),
     surface_(surface),
     coeffsDict_(subDict(type + "Coeffs")),
-    refinementFactor_(readScalar(lookup("refinementFactor")))
+    refinementFactor_
+    (
+        lookupOrDefault<scalar>("refinementFactor", 1.0)
+    )
 {}
 
 
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/uniformValue/uniformValue.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/uniformValue/uniformValue.H
index 9e54e86b48..3677837933 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/uniformValue/uniformValue.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/uniformValue/uniformValue.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation 
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index 09d27dfa36..de23690ea4 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -44,6 +44,142 @@ const Foam::scalar Foam::conformalVoronoiMesh::tolParallel = 1e-3;
 
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
+Foam::scalar Foam::conformalVoronoiMesh::requiredSize
+(
+    const Foam::point& pt
+) const
+{
+    pointIndexHit surfHit;
+    label hitSurface;
+
+    DynamicList<scalar> cellSizeHits;
+
+    geometryToConformTo_.findSurfaceNearest
+    (
+        pt,
+        sqr(GREAT),
+        surfHit,
+        hitSurface
+    );
+
+    if (!surfHit.hit())
+    {
+        FatalErrorIn
+        (
+            "Foam::tensor Foam::conformalVoronoiMesh::requiredAlignment"
+        )
+            << "findSurfaceNearest did not find a hit across the surfaces."
+            << exit(FatalError) << endl;
+    }
+
+    cellSizeHits.append(cellSizeControl().cellSize(pt));
+
+    // Primary alignment
+
+    vectorField norm(1);
+
+    allGeometry_[hitSurface].getNormal
+    (
+        List<pointIndexHit>(1, surfHit),
+        norm
+    );
+
+    const vector np = norm[0];
+
+    // Generate equally spaced 'spokes' in a circle normal to the
+    // direction from the vertex to the closest point on the surface
+    // and look for a secondary intersection.
+
+    const vector d = surfHit.hitPoint() - pt;
+
+    const tensor Rp = rotationTensor(vector(0,0,1), np);
+
+    const label s = cvMeshControls().alignmentSearchSpokes();
+
+    const scalar spanMag = geometryToConformTo_.globalBounds().mag();
+
+    scalar totalDist = 0;
+
+    for (label i = 0; i < s; i++)
+    {
+        vector spoke
+        (
+            Foam::cos(i*constant::mathematical::twoPi/s),
+            Foam::sin(i*constant::mathematical::twoPi/s),
+            0
+        );
+
+        spoke *= spanMag;
+
+        spoke = Rp & spoke;
+
+        pointIndexHit spokeHit;
+
+        label spokeSurface = -1;
+
+        // internal spoke
+
+        geometryToConformTo_.findSurfaceNearestIntersection
+        (
+            pt,
+            pt + spoke,
+            spokeHit,
+            spokeSurface
+        );
+
+        if (spokeHit.hit())
+        {
+            const Foam::point& hitPt = spokeHit.hitPoint();
+
+            scalar spokeHitDistance = mag(hitPt - pt);
+
+            cellSizeHits.append
+            (
+                cellSizeControl().cellSize(hitPt)
+            );
+
+            totalDist += spokeHitDistance;
+        }
+
+        //external spoke
+
+        Foam::point mirrorPt = pt + 2*d;
+
+        geometryToConformTo_.findSurfaceNearestIntersection
+        (
+            mirrorPt,
+            mirrorPt + spoke,
+            spokeHit,
+            spokeSurface
+        );
+
+        if (spokeHit.hit())
+        {
+            const Foam::point& hitPt = spokeHit.hitPoint();
+
+            scalar spokeHitDistance = mag(hitPt - mirrorPt);
+
+            cellSizeHits.append
+            (
+                cellSizeControl().cellSize(hitPt)
+            );
+
+            totalDist += spokeHitDistance;
+        }
+    }
+
+    scalar cellSize = 0;
+
+    forAll(cellSizeHits, hitI)
+    {
+        cellSize += cellSizeHits[hitI];
+    }
+
+    return cellSize/cellSizeHits.size();
+    //return cellSizeControl().cellSize(pt);
+}
+
+
 Foam::tensor Foam::conformalVoronoiMesh::requiredAlignment
 (
     const Foam::point& pt
@@ -828,10 +964,7 @@ void Foam::conformalVoronoiMesh::storeSizesAndAlignments
     {
         sizeAndAlignmentLocations_[i] = topoint(*pit);
 
-        storedSizes_[i] = cellSizeControl().cellSize
-        (
-            sizeAndAlignmentLocations_[i]
-        );
+        storedSizes_[i] = requiredSize(sizeAndAlignmentLocations_[i]);
 
         storedAlignments_[i] = requiredAlignment(sizeAndAlignmentLocations_[i]);
 
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
index 81ce496c71..f9648dc1e7 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
@@ -268,6 +268,9 @@ private:
         //- Return the local maximum surface protrusion distance
         inline scalar maxSurfaceProtrusion(const Foam::point& pt) const;
 
+        //- Return the required cell size at the given location
+        scalar requiredSize(const Foam::point& pt) const;
+
         //- Return the required alignment directions at the given location
         tensor requiredAlignment(const Foam::point& pt) const;
 
@@ -533,6 +536,11 @@ private:
             const Delaunay::Finite_edges_iterator& eit
         ) const;
 
+        boolList dualFaceBoundaryPoints
+        (
+            const Delaunay::Finite_edges_iterator& eit
+        ) const;
+
         //- Finds the maximum filterCount of the dual vertices
         //  (Delaunay cells) that form the dual face produced by the
         //  supplied edge
@@ -839,6 +847,13 @@ private:
             const Delaunay::Finite_facets_iterator& fit
         ) const;
 
+        //- Merge adjacent edges that are not attached to other faces
+        label mergeNearlyParallelEdges
+        (
+            const pointField& pts,
+            const scalar maxCosAngle
+        );
+
         //- Merge vertices that are very close together
         void mergeCloseDualVertices
         (
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
index 5586446cb8..e0fe697e88 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
@@ -212,6 +212,20 @@ void Foam::conformalVoronoiMesh::calcDualMesh
                     Info<< nl << "Collapsing unnecessary faces" << endl;
 
                     collapseFaces(points, boundaryPts, deferredCollapseFaces);
+
+                    const scalar maxCosAngle
+                        = cos(degToRad(cvMeshControls().edgeMergeAngle()));
+
+                    Info<< nl << "Merging adjacent edges which have an angle "
+                        << "of greater than "
+                        << cvMeshControls().edgeMergeAngle() << ": " << endl;
+
+                    label nRemovedEdges =
+                        mergeNearlyParallelEdges(points, maxCosAngle);
+
+                    reduce(nRemovedEdges, sumOp<label>());
+
+                    Info<< "    Merged " << nRemovedEdges << " edges" << endl;
                 }
 
                 labelHashSet wrongFaces = checkPolyMeshQuality(points);
@@ -652,6 +666,113 @@ Foam::label Foam::conformalVoronoiMesh::mergeCloseDualVertices
 }
 
 
+Foam::label Foam::conformalVoronoiMesh::mergeNearlyParallelEdges
+(
+    const pointField& pts,
+    const scalar maxCosAngle
+)
+{
+    List<HashSet<label> > pointFaceCount(number_of_cells());
+    labelList pointNeighbour(number_of_cells(), -1);
+
+    Map<label> dualPtIndexMap;
+
+    for
+    (
+        Delaunay::Finite_edges_iterator eit = finite_edges_begin();
+        eit != finite_edges_end();
+        ++eit
+    )
+    {
+        Cell_handle c = eit->first;
+        Vertex_handle vA = c->vertex(eit->second);
+        Vertex_handle vB = c->vertex(eit->third);
+
+        if (isBoundaryDualFace(eit))
+        {
+            const face f = buildDualFace(eit);
+
+            forAll(f, pI)
+            {
+                const label pIndex = f[pI];
+
+                const label prevPointI = f.prevLabel(pI);
+                const label nextPointI = f.nextLabel(pI);
+
+                pointFaceCount[pIndex].insert(prevPointI);
+                pointFaceCount[pIndex].insert(nextPointI);
+                pointNeighbour[pIndex] = nextPointI;
+            }
+        }
+        else if
+        (
+            vA->internalOrBoundaryPoint()
+         || vB->internalOrBoundaryPoint()
+        )
+        {
+            const face f = buildDualFace(eit);
+            const boolList faceBoundaryPoints = dualFaceBoundaryPoints(eit);
+
+            forAll(f, pI)
+            {
+                const label pIndex = f[pI];
+
+                const label prevPointI = f.prevLabel(pI);
+                const label nextPointI = f.nextLabel(pI);
+
+                pointFaceCount[pIndex].insert(prevPointI);
+                pointFaceCount[pIndex].insert(nextPointI);
+
+                if (faceBoundaryPoints[pI] == false)
+                {
+                    pointNeighbour[pIndex] = nextPointI;
+                }
+                else
+                {
+                    if (faceBoundaryPoints[prevPointI] == true)
+                    {
+                        pointNeighbour[pIndex] = prevPointI;
+                    }
+                    else if (faceBoundaryPoints[nextPointI] == true)
+                    {
+                        pointNeighbour[pIndex] = nextPointI;
+                    }
+                    else
+                    {
+                        pointNeighbour[pIndex] = pIndex;
+                    }
+                }
+            }
+        }
+    }
+
+    forAll(pointFaceCount, pI)
+    {
+        if (pointFaceCount[pI].size() == 2)
+        {
+            List<vector> edges(2, vector(0, 0, 0));
+
+            label count = 0;
+            forAllConstIter(HashSet<label>, pointFaceCount[pI], iter)
+            {
+                edges[count] = pts[pI] - pts[iter.key()];
+                edges[count] /= mag(edges[count]) + VSMALL;
+                count++;
+            }
+
+            if (mag(edges[0] & edges[1]) > maxCosAngle)
+            {
+                dualPtIndexMap.insert(pI, pointNeighbour[pI]);
+            }
+        }
+    }
+
+    reindexDualVertices(dualPtIndexMap);
+
+    return dualPtIndexMap.size();
+}
+
+
 void Foam::conformalVoronoiMesh::smoothSurface
 (
     pointField& pts,
@@ -696,7 +817,6 @@ void Foam::conformalVoronoiMesh::smoothSurface
     } while (nCollapsedFaces > 0);
 
     // Force all points of boundary faces to be on the surface
-
     for
     (
         Delaunay::Finite_cells_iterator cit = finite_cells_begin();
@@ -917,6 +1037,7 @@ void Foam::conformalVoronoiMesh::collapseFaces
 
         mergeCloseDualVertices(pts, boundaryPts);
 
+
         if (nCollapsedFaces > 0)
         {
             Info<< "    Collapsed " << nCollapsedFaces << " faces" << endl;
@@ -932,6 +1053,7 @@ void Foam::conformalVoronoiMesh::collapseFaces
         }
 
     } while (nCollapsedFaces > 0);
+
 }
 
 
@@ -1064,6 +1186,25 @@ Foam::conformalVoronoiMesh::collapseFace
 
     bool allowEarlyCollapseToPoint = true;
 
+
+//    // Quick exit
+//    label smallEdges = 0;
+//    const edgeList& fEdges = f.edges();
+//    forAll(fEdges, eI)
+//    {
+//        const edge& e = fEdges[eI];
+//
+//        if (e.mag(pts) < 0.2*targetFaceSize)
+//        {
+//            smallEdges++;
+//        }
+//    }
+//    if (smallEdges == 0)
+//    {
+//        return fcmNone;
+//    }
+
+
     // if (maxFC > cvMeshControls().filterCountSkipThreshold() - 3)
     // {
     //     limitToQuadsOrTris = true;
@@ -1071,6 +1212,7 @@ Foam::conformalVoronoiMesh::collapseFace
     //     allowEarlyCollapseToPoint = false;
     // }
 
+
     collapseSizeLimitCoeff *= pow
     (
         cvMeshControls().filterErrorReductionCoeff(),
@@ -1158,6 +1300,91 @@ Foam::conformalVoronoiMesh::collapseFace
         }
     }
 
+
+    scalar maxDist = 0;
+    scalar minDist = GREAT;
+//
+//    if (f.size() <= 3)
+//    {
+//        const edgeList& fEdges = f.edges();
+//
+//        forAll(fEdges, eI)
+//        {
+//            const edge& e = fEdges[eI];
+//            const scalar d = e.mag(pts);
+//
+//            if (d > maxDist)
+//            {
+//                maxDist = d;
+//                collapseAxis = e.vec(pts);
+//            }
+//            else if (d < minDist && d != 0)
+//            {
+//                minDist = d;
+//            }
+//        }
+//    }
+//    else
+//    {
+//        forAll(f, pI)
+//        {
+//            for (label i = pI + 1; i < f.size(); ++i)
+//            {
+//                if
+//                (
+//                    f[i] != f.nextLabel(pI)
+//                 && f[i] != f.prevLabel(pI)
+//                )
+//                {
+//                    scalar d = mag(pts[f[pI]] - pts[f[i]]);
+//
+//                    if (d > maxDist)
+//                    {
+//                        maxDist = d;
+//                        collapseAxis = pts[f[pI]] - pts[f[i]];
+//                    }
+//                    else if (d < minDist && d != 0)
+//                    {
+//                        minDist = d;
+//                    }
+//                }
+//            }
+//        }
+//    }
+//
+//    const edgeList& fEdges = f.edges();
+//
+//    scalar perimeter = 0;
+//
+//    forAll(fEdges, eI)
+//    {
+//        const edge& e = fEdges[eI];
+//        const scalar d = e.mag(pts);
+//
+//        perimeter += d;
+//
+////        collapseAxis += e.vec(pts);
+//
+//        if (d > maxDist)
+//        {
+//            collapseAxis = e.vec(pts);
+//            maxDist = d;
+//        }
+//        else if (d < minDist && d != 0)
+//        {
+//            minDist = d;
+//        }
+//    }
+//
+//    collapseAxis /= mag(collapseAxis);
+//
+////    Info<< f.size() << " " << minDist << " " << maxDist << " | "
+////        << collapseAxis << endl;
+//
+//    aspectRatio = maxDist/minDist;
+//
+//    aspectRatio = min(aspectRatio, sqr(perimeter)/(16.0*fA));
+
     if (magSqr(collapseAxis) < VSMALL)
     {
         WarningIn
@@ -1170,9 +1397,9 @@ Foam::conformalVoronoiMesh::collapseFace
         // Output face and collapse axis for visualisation
 
         Pout<< "# Aspect ratio = " << aspectRatio << nl
-            << "# inertia = " << J << nl
-            << "# determinant = " << detJ << nl
-            << "# eigenvalues = " << eigenValues(J) << nl
+//            << "# inertia = " << J << nl
+//            << "# determinant = " << detJ << nl
+//            << "# eigenvalues = " << eigenValues(J) << nl
             << "# collapseAxis = " << collapseAxis << nl
             << "# facePts = " << facePts << nl
             << endl;
@@ -1280,8 +1507,6 @@ Foam::conformalVoronoiMesh::collapseFace
     {
         scalar guardFraction = cvMeshControls().edgeCollapseGuardFraction();
 
-        cvMeshControls().maxCollapseFaceToPointSideLengthCoeff();
-
         if
         (
             allowEarlyCollapseToPoint
@@ -1337,6 +1562,8 @@ Foam::conformalVoronoiMesh::collapseFace
             Foam::point collapseToPt =
                 collapseAxis*(sum(dNeg)/dNeg.size() - dShift) + fC;
 
+//            DynamicList<label> faceBoundaryPts(f.size());
+
             forAll(facePtsNeg, fPtI)
             {
                 if (boundaryPts[facePtsNeg[fPtI]] == true)
@@ -1351,9 +1578,36 @@ Foam::conformalVoronoiMesh::collapseFace
                     collapseToPt = pts[collapseToPtI];
 
                     break;
+
+//                    faceBoundaryPts.append(facePtsNeg[fPtI]);
                 }
             }
 
+//            if (!faceBoundaryPts.empty())
+//            {
+//                if (faceBoundaryPts.size() == 2)
+//                {
+//                    collapseToPtI = faceBoundaryPts[0];
+//
+//                    collapseToPt =
+//                        0.5
+//                       *(
+//                            pts[faceBoundaryPts[0]]
+//                          + pts[faceBoundaryPts[1]]
+//                        );
+//                }
+//                else if (faceBoundaryPts.size() < f.size())
+//                {
+//                    face bFace(faceBoundaryPts);
+//
+//                    collapseToPtI = faceBoundaryPts.first();
+//
+//                    collapseToPt = bFace.centre(pts);
+//                }
+//            }
+//
+//            faceBoundaryPts.clear();
+
             // ...otherwise arbitrarily choosing the most distant
             // point as the index to collapse to.
 
@@ -1384,9 +1638,34 @@ Foam::conformalVoronoiMesh::collapseFace
                     collapseToPt = pts[collapseToPtI];
 
                     break;
+
+//                    faceBoundaryPts.append(facePtsNeg[fPtI]);
                 }
             }
 
+//            if (!faceBoundaryPts.empty())
+//            {
+//                if (faceBoundaryPts.size() == 2)
+//                {
+//                    collapseToPtI = faceBoundaryPts[0];
+//
+//                    collapseToPt =
+//                        0.5
+//                       *(
+//                            pts[faceBoundaryPts[0]]
+//                          + pts[faceBoundaryPts[1]]
+//                        );
+//                }
+//                else if (faceBoundaryPts.size() < f.size())
+//                {
+//                    face bFace(faceBoundaryPts);
+//
+//                    collapseToPtI = faceBoundaryPts.first();
+//
+//                    collapseToPt = bFace.centre(pts);
+//                }
+//            }
+
             // ...otherwise arbitrarily choosing the most distant
             // point as the index to collapse to.
 
@@ -1406,6 +1685,8 @@ Foam::conformalVoronoiMesh::collapseFace
 
             Foam::point collapseToPt = fC;
 
+            DynamicList<label> faceBoundaryPts(f.size());
+
             forAll(facePts, fPtI)
             {
                 if (boundaryPts[facePts[fPtI]] == true)
@@ -1415,11 +1696,32 @@ Foam::conformalVoronoiMesh::collapseFace
                     // use the first boundary point encountered if
                     // there are multiple boundary points.
 
-                    collapseToPtI = facePts[fPtI];
+//                    collapseToPtI = facePts[fPtI];
+//
+//                    collapseToPt = pts[collapseToPtI];
+//
+//                    break;
 
-                    collapseToPt = pts[collapseToPtI];
+                    faceBoundaryPts.append(facePts[fPtI]);
+                }
+            }
 
-                    break;
+            if (!faceBoundaryPts.empty())
+            {
+                if (faceBoundaryPts.size() == 2)
+                {
+                    collapseToPtI = faceBoundaryPts[0];
+
+                    collapseToPt =
+                        0.5*(pts[faceBoundaryPts[0]] + pts[faceBoundaryPts[1]]);
+                }
+                else if (faceBoundaryPts.size() < f.size())
+                {
+                    face bFace(faceBoundaryPts);
+
+                    collapseToPtI = faceBoundaryPts.first();
+
+                    collapseToPt = bFace.centre(pts);
                 }
             }
 
@@ -1791,9 +2093,9 @@ void Foam::conformalVoronoiMesh::checkCellSizing()
         }
 
         Info<< "    Automatically re-sizing " << cellsToResize.size()
-            << " cells that are attached to the bad faces: DISABLED" << endl;
+            << " cells that are attached to the bad faces: " << endl;
 
-        //cellSizeControl_.setCellSizes(cellsToResize);
+        cellSizeControl_.setCellSizes(cellsToResize);
     }
 
     timeCheck("End of Cell Sizing");
@@ -1858,7 +2160,10 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::findOffsetPatchFaces
         }
     }
 
-    offsetBoundaryCells.write();
+    if (cvMeshControls().objOutput())
+    {
+        offsetBoundaryCells.write();
+    }
 
     return offsetBoundaryCells;
 }
@@ -2092,13 +2397,7 @@ void Foam::conformalVoronoiMesh::indexDualVertices
 
             pts[dualVertI] = cit->dual();
 
-            if
-            (
-                !cit->vertex(0)->internalOrBoundaryPoint()
-             || !cit->vertex(1)->internalOrBoundaryPoint()
-             || !cit->vertex(2)->internalOrBoundaryPoint()
-             || !cit->vertex(3)->internalOrBoundaryPoint()
-            )
+            if (cit->boundaryDualVertex())
             {
                 // This is a boundary dual vertex
                 boundaryPts[dualVertI] = true;
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
index bee04b2fbd..bf077cece0 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
@@ -440,6 +440,10 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
                         existingSurfacePtLocations
                     );
                 }
+//                else
+//                {
+//                    vit->setInternal();
+//                }
             }
             else if (vit->ppSlave() || vit->referredExternal())
             {
@@ -477,6 +481,10 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
                         existingSurfacePtLocations
                     );
                 }
+//                else
+//                {
+//                    vit->setInternal();
+//                }
             }
         }
 
@@ -719,12 +727,15 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections
 
                         const Foam::point& p = infoList[hitI].hitPoint();
 
-                        const scalar separationDistance
-                            = mag(p - info.hitPoint());
+                        const scalar separationDistance =
+                            mag(p - info.hitPoint());
 
-                        const scalar minSepDist
-                            = sqr(cvMeshControls().removalDistCoeff()
-                             *targetCellSize(p));
+                        const scalar minSepDist =
+                            sqr
+                            (
+                                cvMeshControls().removalDistCoeff()
+                               *targetCellSize(p)
+                            );
 
                         // Reject the point if it is too close to another
                         // surface point.
@@ -2694,7 +2705,7 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
                             // NEED TO REMOVE FROM THE SURFACE TREE...
                             surfacePtLocationTreePtr_().remove
                             (
-                                existingSurfacePtLocations.size()
+                                existingSurfacePtLocations.size() - 1
                             );
                         }
 
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
index 1e24b7d8ef..43d8c077fc 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
@@ -366,6 +366,49 @@ inline bool Foam::conformalVoronoiMesh::isBoundaryDualFace
 }
 
 
+inline Foam::List<bool> Foam::conformalVoronoiMesh::dualFaceBoundaryPoints
+(
+    const Delaunay::Finite_edges_iterator& eit
+) const
+{
+    Cell_circulator ccStart = incident_cells(*eit);
+    Cell_circulator cc1 = ccStart;
+    Cell_circulator cc2 = cc1;
+
+    // Advance the second circulator so that it always stays on the next
+    // cell around the edge;
+    cc2++;
+
+    DynamicList<bool> tmpFaceBoundaryPoints;
+
+    do
+    {
+        label cc1I = cc1->cellIndex();
+
+        label cc2I = cc2->cellIndex();
+
+        if (cc1I != cc2I)
+        {
+            if (cc1->boundaryDualVertex())
+            {
+                tmpFaceBoundaryPoints.append(true);
+            }
+            else
+            {
+                tmpFaceBoundaryPoints.append(false);
+            }
+        }
+
+        cc1++;
+
+        cc2++;
+
+    } while (cc1 != ccStart);
+
+    return tmpFaceBoundaryPoints;
+}
+
+
 inline Foam::List<Foam::label> Foam::conformalVoronoiMesh::processorsAttached
 (
     const Delaunay::Finite_facets_iterator& fit
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.C
index c99f9324ff..996aa71abf 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.C
@@ -291,6 +291,11 @@ Foam::cvControls::cvControls
         filteringDict.lookup("mergeClosenessCoeff")
     );
 
+    edgeMergeAngle_ = readScalar
+    (
+        filteringDict.lookup("edgeMergeAngle")
+    );
+
     continueFilteringOnBadInitialPolyMesh_ = Switch
     (
         filteringDict.lookupOrDefault<Switch>
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.H
index dae83ed06d..58dd06fc74 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.H
@@ -234,6 +234,11 @@ class cvControls
         //  being merged, fraction of the local target cell size
         scalar mergeClosenessCoeff_;
 
+        //- If the angle between two dual edges that are connected by a single
+        //  point is less than this angle, then the edges will be merged into a
+        //  single edge.
+        scalar edgeMergeAngle_;
+
         //- If the mesh quality criteria cannot be satisfied, continue
         //  with filtering anyway?
         Switch continueFilteringOnBadInitialPolyMesh_;
@@ -405,6 +410,9 @@ public:
             //- Return the mergeClosenessCoeff
             inline scalar mergeClosenessCoeff() const;
 
+            //- Return the edgeMergeAngle
+            inline scalar edgeMergeAngle() const;
+
             //- Return the continueFilteringOnBadInitialPolyMesh Switch
             inline Switch continueFilteringOnBadInitialPolyMesh() const;
 
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControlsI.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControlsI.H
index ce45e5f522..50df7e8f58 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControlsI.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControlsI.H
@@ -166,6 +166,12 @@ inline Foam::scalar Foam::cvControls::mergeClosenessCoeff() const
 }
 
 
+inline Foam::scalar Foam::cvControls::edgeMergeAngle() const
+{
+    return edgeMergeAngle_;
+}
+
+
 inline Foam::Switch
 Foam::cvControls::continueFilteringOnBadInitialPolyMesh() const
 {
diff --git a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C
index a0873ab4c9..b1bcf33d82 100644
--- a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C
+++ b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C
@@ -36,7 +36,6 @@ Description
       (nonuniformTransform)cyclic \<zoneA\>_\<zoneB\>
     - extrudes into master direction (i.e. away from the owner cell
       if flipMap is false)
-    - not parallel
 
 \verbatim
 
diff --git a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/extrude2DMesh/extrude2DMesh.C b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/extrude2DMesh/extrude2DMesh.C
index f7c1026542..0bbe7d926a 100644
--- a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/extrude2DMesh/extrude2DMesh.C
+++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/extrude2DMesh/extrude2DMesh.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/extrude2DMesh/extrude2DMesh.H b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/extrude2DMesh/extrude2DMesh.H
index 0733b5e236..3b6efd267e 100644
--- a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/extrude2DMesh/extrude2DMesh.H
+++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/extrude2DMesh/extrude2DMesh.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.C b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.C
index 30c7a0334a..7979fa966a 100644
--- a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.C
+++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -330,4 +330,4 @@ void Foam::patchToPoly2DMesh::createMesh()
 // * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
 
 
-// ************************************************************************* // 
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.H b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.H
index a07f73f0d4..ea27dc6f53 100644
--- a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.H
+++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -162,4 +162,4 @@ public:
 
 #endif
 
-// ************************************************************************* // 
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMeshApp.C b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMeshApp.C
index 4bb6b97b16..5381c0b57a 100644
--- a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMeshApp.C
+++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMeshApp.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
index a01fa988d2..28930c573f 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
@@ -82,7 +82,7 @@ Foam::label Foam::checkTopology
                 if (cFaces[i] < 0 || cFaces[i] >= mesh.nFaces())
                 {
                     cells.insert(cellI);
-                    break;  
+                    break;
                 }
             }
         }
diff --git a/applications/utilities/mesh/manipulation/createBaffles/createBaffles.C b/applications/utilities/mesh/manipulation/createBaffles/createBaffles.C
index b3d9d3f515..6fb57aa027 100644
--- a/applications/utilities/mesh/manipulation/createBaffles/createBaffles.C
+++ b/applications/utilities/mesh/manipulation/createBaffles/createBaffles.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -142,14 +142,6 @@ int main(int argc, char *argv[])
         "additionalPatches",
         "((master2 slave2) .. (masterN slaveN))"
     );
-    argList::addBoolOption("internalFacesOnly");
-
-    argList::addOption
-    (
-        "additionalPatches",
-        "(patch2 .. patchN)",
-        "specify additional patches for creating baffles"
-    );
     argList::addBoolOption
     (
         "internalFacesOnly",
diff --git a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
index b34676fac5..d129efde6d 100644
--- a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
+++ b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -2194,6 +2194,8 @@ int main(int argc, char *argv[])
 
             label regionI = -1;
 
+            (void)mesh.tetBasePtIs();
+
             label cellI = mesh.findCell(insidePoint);
 
             Info<< nl << "Found point " << insidePoint << " in cell " << cellI
diff --git a/applications/utilities/postProcessing/sampling/sample/sampleDict b/applications/utilities/postProcessing/sampling/sample/sampleDict
index b55b190bf7..7a9b6f7f0a 100644
--- a/applications/utilities/postProcessing/sampling/sample/sampleDict
+++ b/applications/utilities/postProcessing/sampling/sample/sampleDict
@@ -47,12 +47,14 @@ formatOptions
 }
 
 // interpolationScheme. choice of
-//      cell          : use cell-centre value only; constant over cells (default)
+//      cell          : use cell-centre value only; constant over cells
+//                      (default)
 //      cellPoint     : use cell-centre and vertex values
 //      cellPointFace : use cell-centre, vertex and face values.
 //      pointMVC      : use point values only (Mean Value Coordinates)
-//      cellPatchConstrained : use cell-centre except on boundary faces where
-//        it uses the boundary value. For use with e.g. patchCloudSet.
+//      cellPatchConstrained : like 'cell' but uses cell-centre except on
+//                             boundary faces where it uses the boundary value.
+//                             For use with e.g. patchCloudSet.
 // 1] vertex values determined from neighbouring cell-centre values
 // 2] face values determined using the current face interpolation scheme
 //    for the field (linear, gamma, etc.)
diff --git a/applications/utilities/surface/surfaceFeatureExtract/Allwmake b/applications/utilities/surface/surfaceFeatureExtract/Allwmake
index e54f15ae28..369428e171 100755
--- a/applications/utilities/surface/surfaceFeatureExtract/Allwmake
+++ b/applications/utilities/surface/surfaceFeatureExtract/Allwmake
@@ -3,30 +3,25 @@ cd ${0%/*} || exit 1    # run from this directory
 
 set -x
 
-if [ ! -e "Make/files" ] || [ ! -e "Make/options" ]
+if [ -n "$CGAL_ARCH_PATH" ]
 then
-    mkdir -p Make
+    echo
+    echo "Compiling surfaceFeatureExtract with CGAL curvature support"
+    echo
 
-    if [ -n "$CGAL_ARCH_PATH" ]
-    then
-        cp -r MakeWithCGAL/* Make
-
-        echo
-        echo Compiling surfaceFeatureExtract with CGAL support for curvature
-        echo
-
-        wmake
-    else
-        cp -r MakeWithoutCGAL/* Make
-
-        echo
-        echo Compiling surfaceFeatureExtract without CGAL support for curvature
-        echo
-
-        wmake
-    fi
+    wmake "ENABLE_CURVATURE=-DENABLE_CURVATURE \
+           EXE_FROUNDING_MATH=-frounding-math \
+           USE_F2C=-DCGAL_USE_F2C \
+           CGAL_LIBDIR=-L$CGAL_ARCH_PATH/lib \
+           LAPACK_LIB=-llapack \
+           BLAS_LIB=-lblas \
+           CGAL_LIB=-lCGAL"
 else
-    echo surfaceFeatureExtract already has a Make folder
+    echo
+    echo "Compiling surfaceFeatureExtract without CGAL curvature support"
+    echo
+
+    wmake
 fi
 
 
diff --git a/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/buildCGALPolyhedron.H b/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/buildCGALPolyhedron.H
index 3001a247d3..30627c7463 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/buildCGALPolyhedron.H
+++ b/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/buildCGALPolyhedron.H
@@ -53,7 +53,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                     Class buildCGALPolyhedron Declaration
+                    Class buildCGALPolyhedron Declaration
 \*---------------------------------------------------------------------------*/
 
 class buildCGALPolyhedron
@@ -80,18 +80,56 @@ public:
     // Constructors
 
         //- Construct with reference to triSurface
-        buildCGALPolyhedron(const triSurface& surf);
+        explicit buildCGALPolyhedron(const triSurface& surf)
+        :
+            CGAL::Modifier_base<HalfedgeDS>(),
+            surf_(surf)
+        {}
 
 
     //- Destructor
-    ~buildCGALPolyhedron();
+    ~buildCGALPolyhedron(){}
 
 
     // Member Operators
 
         //- operator() of this `modifier' called by delegate function of
         //  Polyhedron
-        void operator()(HalfedgeDS& hds);
+        void operator()(HalfedgeDS& hds)
+        {
+            typedef HalfedgeDS::Traits     Traits;
+            typedef Traits::Point_3        Point;
+
+            // Postcondition: `hds' is a valid polyhedral surface.
+            CGAL::Polyhedron_incremental_builder_3<HalfedgeDS> B(hds, false);
+
+            B.begin_surface
+            (
+                surf_.points().size(),      // n points
+                surf_.size(),               // n facets
+                2*surf_.edges().size()      // n halfedges
+            );
+
+            forAll(surf_.points(), pI)
+            {
+                const Foam::point& p = surf_.points()[pI];
+
+                B.add_vertex(Point(p.x(), p.y(), p.z()));
+            }
+
+            forAll(surf_, fI)
+            {
+                B.begin_facet();
+
+                B.add_vertex_to_facet(surf_[fI][0]);
+                B.add_vertex_to_facet(surf_[fI][1]);
+                B.add_vertex_to_facet(surf_[fI][2]);
+
+                B.end_facet();
+            }
+
+            B.end_surface();
+        }
 };
 
 
diff --git a/applications/utilities/surface/surfaceFeatureExtract/Make/options b/applications/utilities/surface/surfaceFeatureExtract/Make/options
index b733d5fde9..b2e5dcf5f7 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/Make/options
+++ b/applications/utilities/surface/surfaceFeatureExtract/Make/options
@@ -1,4 +1,11 @@
+include $(GENERAL_RULES)/CGAL
+
 EXE_INC = \
+    ${ENABLE_CURVATURE}\
+    ${EXE_FROUNDING_MATH} \
+    ${USE_F2C} \
+    ${CGAL_INC} \
+    -ICGALPolyhedron \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/edgeMesh/lnInclude \
@@ -7,6 +14,11 @@ EXE_INC = \
     -I$(LIB_SRC)/sampling/lnInclude
 
 EXE_LIBS = \
+    $(CGAL_LIBS) \
+    ${CGAL_LIBDIR} \
+    ${LAPACK_LIB} \
+    ${BLAS_LIB} \
+    ${CGAL_LIB} \
     -lmeshTools \
     -ledgeMesh \
     -ltriSurface \
diff --git a/applications/utilities/surface/surfaceFeatureExtract/MakeWithCGAL/files b/applications/utilities/surface/surfaceFeatureExtract/MakeWithCGAL/files
deleted file mode 100644
index 0c0f6f7966..0000000000
--- a/applications/utilities/surface/surfaceFeatureExtract/MakeWithCGAL/files
+++ /dev/null
@@ -1,4 +0,0 @@
-surfaceFeatureExtract.C
-CGALPolyhedron/buildCGALPolyhedron.C
-
-EXE = $(FOAM_APPBIN)/surfaceFeatureExtract
diff --git a/applications/utilities/surface/surfaceFeatureExtract/MakeWithCGAL/options b/applications/utilities/surface/surfaceFeatureExtract/MakeWithCGAL/options
deleted file mode 100644
index 60833572b7..0000000000
--- a/applications/utilities/surface/surfaceFeatureExtract/MakeWithCGAL/options
+++ /dev/null
@@ -1,29 +0,0 @@
-EXE_FROUNDING_MATH = -frounding-math
-EXE_NDEBUG = -DNDEBUG
-USE_F2C = -DCGAL_USE_F2C
-include $(GENERAL_RULES)/CGAL
-
-EXE_INC = \
-    -DENABLE_CURVATURE \
-    ${EXE_FROUNDING_MATH} \
-    ${EXE_NDEBUG} \
-    ${USE_F2C} \
-    ${CGAL_INC} \
-    -ICGALPolyhedron \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/edgeMesh/lnInclude \
-    -I$(LIB_SRC)/triSurface/lnInclude \
-    -I$(LIB_SRC)/surfMesh/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude
-
-EXE_LIBS = \
-    $(CGAL_LIBS) \
-    -L$(CGAL_ARCH_PATH)/lib \
-    -llapack \
-    -lblas \
-    -lCGAL \
-    -lmeshTools \
-    -ledgeMesh \
-    -ltriSurface \
-    -lsampling
diff --git a/applications/utilities/surface/surfaceFeatureExtract/MakeWithoutCGAL/files b/applications/utilities/surface/surfaceFeatureExtract/MakeWithoutCGAL/files
deleted file mode 100644
index a57125cd6a..0000000000
--- a/applications/utilities/surface/surfaceFeatureExtract/MakeWithoutCGAL/files
+++ /dev/null
@@ -1,3 +0,0 @@
-surfaceFeatureExtract.C
-
-EXE = $(FOAM_APPBIN)/surfaceFeatureExtract
diff --git a/applications/utilities/surface/surfaceFeatureExtract/MakeWithoutCGAL/options b/applications/utilities/surface/surfaceFeatureExtract/MakeWithoutCGAL/options
deleted file mode 100644
index b733d5fde9..0000000000
--- a/applications/utilities/surface/surfaceFeatureExtract/MakeWithoutCGAL/options
+++ /dev/null
@@ -1,13 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/edgeMesh/lnInclude \
-    -I$(LIB_SRC)/triSurface/lnInclude \
-    -I$(LIB_SRC)/surfMesh/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude
-
-EXE_LIBS = \
-    -lmeshTools \
-    -ledgeMesh \
-    -ltriSurface \
-    -lsampling
diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
index c2baec519e..5b2e78ca96 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -74,30 +74,6 @@ scalarField calcCurvature(const triSurface& surf)
 
     // The rest of this function adapted from
     //     CGAL-3.7/examples/Jet_fitting_3/Mesh_estimation.cpp
-    // Licensed under CGAL-3.7/LICENSE.FREE_USE
-
-    // Copyright (c) 1996,1997,1998,1999,2000,2001,2002,2003,2004,2005,2006,2007
-    // Utrecht University (The Netherlands), ETH Zurich (Switzerland), Freie
-    // Universitaet Berlin (Germany), INRIA Sophia-Antipolis (France),
-    // Martin-Luther-University Halle-Wittenberg (Germany), Max-Planck-Institute
-    // Saarbruecken (Germany), RISC Linz (Austria), and Tel-Aviv University
-    // (Israel).  All rights reserved.
-
-    // Permission is hereby granted, free of charge, to any person obtaining a
-    // copy of this software and associated documentation files (the
-    // "Software"), to deal in the Software without restriction, including
-    // without limitation the rights to use, copy, modify, merge, publish,
-    // distribute, sublicense, and/or sell copies of the Software, and to permit
-    // persons to whom the Software is furnished to do so, subject to the
-    // following conditions:
-
-    // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
-    // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
-    // MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
-    // IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
-    // CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT
-    // OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR
-    // THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 
      //Vertex property map, with std::map
     typedef std::map<Vertex*, int> Vertex2int_map_type;
@@ -181,12 +157,12 @@ scalarField calcCurvature(const triSurface& surf)
 //                 << monge_fit.condition_number() << nl << std::endl;
 
         // Use the maximum curvature to give smaller cell sizes later.
-        k[vertI++]
-            = max
-              (
-                   mag(monge_form.principal_curvatures(0)),
-                   mag(monge_form.principal_curvatures(1))
-              );
+        k[vertI++] =
+            max
+            (
+                mag(monge_form.principal_curvatures(0)),
+                mag(monge_form.principal_curvatures(1))
+            );
     }
 
     return k;
@@ -198,8 +174,10 @@ bool edgesConnected(const edge& e1, const edge& e2)
 {
     if
     (
-        e1.start() == e2.start() || e1.start() == e2.end()
-     || e1.end() == e2.start() || e1.end() == e2.end()
+        e1.start() == e2.start()
+     || e1.start() == e2.end()
+     || e1.end() == e2.start()
+     || e1.end() == e2.end()
     )
     {
         return true;
@@ -289,8 +267,8 @@ scalar calcProximityOfFeatureEdges
                     // Don't refine if the edges are connected to each other
                     if (!edgesConnected(e1, e2))
                     {
-                        scalar curDist
-                            = mag(pHit1.hitPoint() - pHit2.hitPoint());
+                        scalar curDist =
+                            mag(pHit1.hitPoint() - pHit2.hitPoint());
 
                         minDist = min(curDist, minDist);
                     }
@@ -503,358 +481,344 @@ int main(int argc, char *argv[])
         "extract and write surface features to file"
     );
     argList::noParallel();
-    argList::validArgs.append("surface");
-    argList::validArgs.append("output set");
 
     argList::addOption
     (
-        "includedAngle",
-        "degrees",
-        "construct feature set from included angle [0..180]"
+        "dict",
+        "word",
+        "specify alternative dictionary for the feature extraction information"
     );
-    argList::addOption
-    (
-        "set",
-        "name",
-        "use existing feature set from file"
-    );
-    argList::addOption
-    (
-        "minLen",
-        "scalar",
-        "remove features shorter than the specified cumulative length"
-    );
-    argList::addOption
-    (
-        "minElem",
-        "int",
-        "remove features with fewer than the specified number of edges"
-    );
-    argList::addOption
-    (
-        "subsetBox",
-        "((x0 y0 z0)(x1 y1 z1))",
-        "remove edges outside specified bounding box"
-    );
-    argList::addOption
-    (
-        "deleteBox",
-        "((x0 y0 z0)(x1 y1 z1))",
-        "remove edges within specified bounding box"
-    );
-    argList::addBoolOption
-    (
-        "writeObj",
-        "write extendedFeatureEdgeMesh obj files"
-    );
-    argList::addBoolOption
-    (
-        "writeVTK",
-        "write extendedFeatureEdgeMesh vtk files"
-    );
-    argList::addBoolOption
-    (
-        "closeness",
-        "span to look for surface closeness"
-    );
-    argList::addOption
-    (
-        "featureProximity",
-        "scalar",
-        "distance to look for close features"
-    );
-    argList::addBoolOption
-    (
-        "writeVTK",
-        "write surface property VTK files"
-    );
-    argList::addBoolOption
-    (
-        "manifoldEdgesOnly",
-        "remove any non-manifold (open or more than two connected faces) edges"
-    );
-    argList::addOption
-    (
-        "plane",
-        "(nx ny nz)(z0 y0 z0)",
-        "use a plane to create feature edges for 2D meshing"
-    );
-
-#   ifdef ENABLE_CURVATURE
-    argList::addBoolOption
-    (
-        "curvature",
-        "calculate curvature and closeness fields"
-    );
-#   endif
-
 
 #   include "setRootCase.H"
 #   include "createTime.H"
 
-    bool writeVTK = args.optionFound("writeVTK");
-
-    bool writeObj = args.optionFound("writeObj");
-
-    bool curvature = args.optionFound("curvature");
-
-    if (curvature && env("FOAM_SIGFPE"))
-    {
-        WarningIn(args.executable())
-            << "Detected floating point exception trapping (FOAM_SIGFPE)."
-            << " This might give" << nl
-            << "    problems when calculating curvature on straight angles"
-            << " (infinite curvature)" << nl
-            << "    Switch it off in case of problems." << endl;
-    }
-
-
-    Info<< "Feature line extraction is only valid on closed manifold surfaces."
-        << endl;
-
-    const fileName surfFileName = args[1];
-    const fileName outFileName  = args[2];
-
-    Info<< "Surface            : " << surfFileName << nl
-        << "Output feature set : " << outFileName << nl
-        << endl;
-
-    fileName sFeatFileName = surfFileName.lessExt().name();
-
-
-    // Read
-    // ~~~~
-
-    triSurface surf(surfFileName);
-
-    Info<< "Statistics:" << endl;
-    surf.writeStats(Info);
-    Info<< endl;
-
-    faceList faces(surf.size());
-
-    forAll(surf, fI)
-    {
-        faces[fI] = surf[fI].triFaceFace();
-    }
-
-    // Either construct features from surface&featureangle or read set.
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    surfaceFeatures set(surf);
-
-    if (args.optionFound("set"))
-    {
-        const fileName setName = args["set"];
-
-        Info<< "Reading existing feature set from file " << setName << endl;
-
-        set = surfaceFeatures(surf, setName);
-    }
-    else if (args.optionFound("includedAngle"))
-    {
-        const scalar includedAngle = args.optionRead<scalar>("includedAngle");
-
-        Info<< "Constructing feature set from included angle " << includedAngle
-            << endl;
-
-        set = surfaceFeatures(surf, includedAngle);
-
-        // Info<< nl << "Writing initial features" << endl;
-        // set.write("initial.fSet");
-        // set.writeObj("initial");
-    }
-    else
-    {
-        FatalErrorIn(args.executable())
-            << "No initial feature set. Provide either one"
-            << " of -set (to read existing set)" << nl
-            << " or -includedAngle (to new set construct from angle)"
-            << exit(FatalError);
-    }
-
-    Info<< nl
-        << "Initial feature set:" << nl
-        << "    feature points : " << set.featurePoints().size() << nl
-        << "    feature edges  : " << set.featureEdges().size() << nl
-        << "    of which" << nl
-        << "        region edges   : " << set.nRegionEdges() << nl
-        << "        external edges : " << set.nExternalEdges() << nl
-        << "        internal edges : " << set.nInternalEdges() << nl
-        << endl;
-
-
-    // Trim set
-    // ~~~~~~~~
-
-    scalar minLen = -GREAT;
-    if (args.optionReadIfPresent("minLen", minLen))
-    {
-        Info<< "Removing features of length < " << minLen << endl;
-    }
-
-    label minElem = 0;
-    if (args.optionReadIfPresent("minElem", minElem))
-    {
-        Info<< "Removing features with number of edges < " << minElem << endl;
-    }
-
-    // Trim away small groups of features
-    if (minElem > 0 || minLen > 0)
-    {
-        set.trimFeatures(minLen, minElem);
-        Info<< endl << "Removed small features" << endl;
-    }
-
-
-    // Subset
-    // ~~~~~~
-
-    // Convert to marked edges, points
-    List<surfaceFeatures::edgeStatus> edgeStat(set.toStatus());
-
-    if (args.optionFound("subsetBox"))
-    {
-        treeBoundBox bb
-        (
-            args.optionLookup("subsetBox")()
-        );
-
-        Info<< "Removing all edges outside bb " << bb << endl;
-        dumpBox(bb, "subsetBox.obj");
-
-        deleteBox
-        (
-            surf,
-            bb,
-            false,
-            edgeStat
-        );
-    }
-    else if (args.optionFound("deleteBox"))
-    {
-        treeBoundBox bb
-        (
-            args.optionLookup("deleteBox")()
-        );
-
-        Info<< "Removing all edges inside bb " << bb << endl;
-        dumpBox(bb, "deleteBox.obj");
-
-        deleteBox
-        (
-            surf,
-            bb,
-            true,
-            edgeStat
-        );
-    }
-
-    if (args.optionFound("manifoldEdgesOnly"))
-    {
-        Info<< "Removing all non-manifold edges" << endl;
-
-        forAll(edgeStat, edgeI)
-        {
-            if (surf.edgeFaces()[edgeI].size() != 2)
-            {
-                edgeStat[edgeI] = surfaceFeatures::NONE;
-            }
-        }
-    }
-
-    if (args.optionFound("plane"))
-    {
-        plane cutPlane(args.optionLookup("plane")());
-
-        deleteEdges
-        (
-            surf,
-            cutPlane,
-            edgeStat
-        );
-
-        Info<< "Only edges that intersect the plane with normal "
-            << cutPlane.normal() << " and base point " << cutPlane.refPoint()
-            << " will be included as feature edges."<< endl;
-    }
-
-
-    surfaceFeatures newSet(surf);
-    newSet.setFromStatus(edgeStat);
-
-    //Info<< endl << "Writing trimmed features to " << outFileName << endl;
-    //newSet.write(outFileName);
-
-    // Info<< endl << "Writing edge objs." << endl;
-    // newSet.writeObj("final");
-
-    Info<< nl
-        << "Final feature set:" << nl
-        << "    feature points : " << newSet.featurePoints().size() << nl
-        << "    feature edges  : " << newSet.featureEdges().size() << nl
-        << "    of which" << nl
-        << "        region edges   : " << newSet.nRegionEdges() << nl
-        << "        external edges : " << newSet.nExternalEdges() << nl
-        << "        internal edges : " << newSet.nInternalEdges() << nl
-        << endl;
-
-    // Extracting and writing a extendedFeatureEdgeMesh
-    extendedFeatureEdgeMesh feMesh
+    word dictName
     (
-        newSet,
-        runTime,
-        sFeatFileName + ".extendedFeatureEdgeMesh"
+        args.optionLookupOrDefault<word>("dict", "surfaceFeatureExtractDict")
     );
 
-    Info<< nl << "Writing extendedFeatureEdgeMesh to " << feMesh.objectPath()
-        << endl;
+    Info<< "Reading " << dictName << nl << endl;
 
-    if (writeObj)
-    {
-        feMesh.writeObj(surfFileName.lessExt().name());
-    }
-
-    feMesh.write();
-
-    // Write a featureEdgeMesh for backwards compatibility
-    {
-        featureEdgeMesh bfeMesh
-        (
-            IOobject
-            (
-                surfFileName.lessExt().name() + ".eMesh",   // name
-                runTime.constant(), // instance
-                "triSurface",
-                runTime,            // registry
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE,
-                false
-            ),
-            feMesh.points(),
-            feMesh.edges()
-        );
-
-        Info<< nl << "Writing featureEdgeMesh to "
-            << bfeMesh.objectPath() << endl;
-
-        bfeMesh.regIOobject::write();
-    }
-
-    triSurfaceMesh searchSurf
+    IOdictionary dict
     (
         IOobject
         (
-            sFeatFileName + ".closeness",
-            runTime.constant(),
-            "triSurface",
+            dictName,
+            runTime.system(),
             runTime,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        surf
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false
+        )
     );
 
+    forAllConstIter(dictionary, dict, iter)
+    {
+        dictionary surfaceDict = dict.subDict(iter().keyword());
+
+        const fileName surfFileName = iter().keyword();
+        const fileName sFeatFileName = surfFileName.lessExt().name();
+
+        Info<< "Surface            : " << surfFileName << nl << endl;
+
+        const Switch writeVTK =
+            surfaceDict.lookupOrAddDefault<Switch>("writeVTK", "off");
+        const Switch writeObj =
+            surfaceDict.lookupOrAddDefault<Switch>("writeObj", "off");
+        const Switch writeFeatureEdgeMesh =
+            surfaceDict.lookupOrAddDefault<Switch>
+            (
+                "writeFeatureEdgeMesh",
+                "off"
+            );
+
+        const Switch curvature =
+            surfaceDict.lookupOrAddDefault<Switch>("curvature", "off");
+        const Switch featureProximity =
+            surfaceDict.lookupOrAddDefault<Switch>("featureProximity", "off");
+        const Switch closeness =
+            surfaceDict.lookupOrAddDefault<Switch>("closeness", "off");
+
+        const word extractionMethod = surfaceDict.lookup("extractionMethod");
+
+
+#ifndef ENABLE_CURVATURE
+        if (curvature)
+        {
+            WarningIn(args.executable())
+                << "Curvature calculation has been requested but "
+                << args.executable() << " has not " << nl
+                << "    been compiled with CGAL. "
+                << "Skipping the curvature calculation." << endl;
+        }
+#else
+        if (curvature && env("FOAM_SIGFPE"))
+        {
+            WarningIn(args.executable())
+                << "Detected floating point exception trapping (FOAM_SIGFPE)."
+                << " This might give" << nl
+                << "    problems when calculating curvature on straight angles"
+                << " (infinite curvature)" << nl
+                << "    Switch it off in case of problems." << endl;
+        }
+#endif
+
+
+        Info<< nl << "Feature line extraction is only valid on closed manifold "
+            << "surfaces." << endl;
+
+        // Read
+        // ~~~~
+
+        triSurface surf("constant/triSurface/" + surfFileName);
+
+        Info<< "Statistics:" << endl;
+        surf.writeStats(Info);
+        Info<< endl;
+
+        faceList faces(surf.size());
+
+        forAll(surf, fI)
+        {
+            faces[fI] = surf[fI].triFaceFace();
+        }
+
+
+        // Either construct features from surface & featureAngle or read set.
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        surfaceFeatures set(surf);
+
+        if (extractionMethod == "extractFromFile")
+        {
+            const fileName featureEdgeFile =
+                surfaceDict.subDict("extractFromFile").lookup
+                (
+                    "featureEdgeFile"
+                );
+
+            edgeMesh eMesh(featureEdgeFile);
+
+            // Sometimes duplicate edges are present. Remove them.
+            eMesh.mergeEdges();
+
+            Info<< nl << "Reading existing feature edges from file "
+                << featureEdgeFile << endl;
+
+            set = surfaceFeatures(surf, eMesh.points(), eMesh.edges());
+        }
+        else if (extractionMethod == "extractFromSurface")
+        {
+            const scalar includedAngle =
+                readScalar
+                (
+                    surfaceDict.subDict("extractFromSurface").lookup
+                    (
+                        "includedAngle"
+                    )
+                );
+
+            Info<< nl << "Constructing feature set from included angle "
+                << includedAngle << endl;
+
+            set = surfaceFeatures(surf, includedAngle);
+        }
+        else
+        {
+            FatalErrorIn(args.executable())
+                << "No initial feature set. Provide either one"
+                << " of extractFromFile (to read existing set)" << nl
+                << " or extractFromSurface (to construct new set from angle)"
+                << exit(FatalError);
+        }
+
+        Info<< nl
+            << "Initial feature set:" << nl
+            << "    feature points : " << set.featurePoints().size() << nl
+            << "    feature edges  : " << set.featureEdges().size() << nl
+            << "    of which" << nl
+            << "        region edges   : " << set.nRegionEdges() << nl
+            << "        external edges : " << set.nExternalEdges() << nl
+            << "        internal edges : " << set.nInternalEdges() << nl
+            << endl;
+
+
+        // Trim set
+        // ~~~~~~~~
+
+        if (surfaceDict.isDict("trimFeatures"))
+        {
+            dictionary trimDict = surfaceDict.subDict("trimFeatures");
+
+            scalar minLen =
+                trimDict.lookupOrAddDefault<scalar>("minLen", -GREAT);
+
+            label minElem = trimDict.lookupOrAddDefault<label>("minElem", 0);
+
+            // Trim away small groups of features
+            if (minElem > 0 || minLen > 0)
+            {
+                Info<< "Removing features of length < "
+                    << minLen << endl;
+                Info<< "Removing features with number of edges < "
+                    << minElem << endl;
+
+                set.trimFeatures(minLen, minElem);
+            }
+        }
+
+
+        // Subset
+        // ~~~~~~
+
+        // Convert to marked edges, points
+        List<surfaceFeatures::edgeStatus> edgeStat(set.toStatus());
+
+        if (surfaceDict.isDict("subsetFeatures"))
+        {
+            dictionary subsetDict = surfaceDict.subDict("subsetFeatures");
+
+            if (subsetDict.found("insideBox"))
+            {
+                treeBoundBox bb(subsetDict.lookup("insideBox")());
+
+                Info<< "Removing all edges outside bb " << bb << endl;
+                dumpBox(bb, "subsetBox.obj");
+
+                deleteBox
+                (
+                    surf,
+                    bb,
+                    false,
+                    edgeStat
+                );
+            }
+            else if (subsetDict.found("outsideBox"))
+            {
+                treeBoundBox bb(subsetDict.lookup("outsideBox")());
+
+                Info<< "Removing all edges inside bb " << bb << endl;
+                dumpBox(bb, "deleteBox.obj");
+
+                deleteBox
+                (
+                    surf,
+                    bb,
+                    true,
+                    edgeStat
+                );
+            }
+
+            const Switch manifoldEdges =
+                subsetDict.lookupOrAddDefault<Switch>("manifoldEdges", "no");
+
+            if (manifoldEdges)
+            {
+                Info<< "Removing all non-manifold edges" << endl;
+
+                forAll(edgeStat, edgeI)
+                {
+                    if (surf.edgeFaces()[edgeI].size() != 2)
+                    {
+                        edgeStat[edgeI] = surfaceFeatures::NONE;
+                    }
+                }
+            }
+
+            if (subsetDict.found("plane"))
+            {
+                plane cutPlane(subsetDict.lookup("plane")());
+
+                deleteEdges
+                (
+                    surf,
+                    cutPlane,
+                    edgeStat
+                );
+
+                Info<< "Only edges that intersect the plane with normal "
+                    << cutPlane.normal()
+                    << " and base point " << cutPlane.refPoint()
+                    << " will be included as feature edges."<< endl;
+            }
+        }
+
+
+        surfaceFeatures newSet(surf);
+        newSet.setFromStatus(edgeStat);
+
+        if (writeObj)
+        {
+            newSet.writeObj("final");
+        }
+
+        Info<< nl
+            << "Final feature set after trimming and subsetting:" << nl
+            << "    feature points : " << newSet.featurePoints().size() << nl
+            << "    feature edges  : " << newSet.featureEdges().size() << nl
+            << "    of which" << nl
+            << "        region edges   : " << newSet.nRegionEdges() << nl
+            << "        external edges : " << newSet.nExternalEdges() << nl
+            << "        internal edges : " << newSet.nInternalEdges() << nl
+            << endl;
+
+        // Extracting and writing a extendedFeatureEdgeMesh
+        extendedFeatureEdgeMesh feMesh
+        (
+            newSet,
+            runTime,
+            sFeatFileName + ".extendedFeatureEdgeMesh"
+        );
+
+        Info<< nl << "Writing extendedFeatureEdgeMesh to "
+            << feMesh.objectPath() << endl;
+
+        if (writeObj)
+        {
+            feMesh.writeObj(surfFileName.lessExt().name());
+        }
+
+        feMesh.write();
+
+        // Write a featureEdgeMesh for backwards compatibility
+        if (writeFeatureEdgeMesh)
+        {
+            featureEdgeMesh bfeMesh
+            (
+                IOobject
+                (
+                    surfFileName.lessExt().name() + ".eMesh",   // name
+                    runTime.constant(),                         // instance
+                    "triSurface",
+                    runTime,                                    // registry
+                    IOobject::NO_READ,
+                    IOobject::AUTO_WRITE,
+                    false
+                ),
+                feMesh.points(),
+                feMesh.edges()
+            );
+
+            Info<< nl << "Writing featureEdgeMesh to "
+                << bfeMesh.objectPath() << endl;
+
+            bfeMesh.regIOobject::write();
+        }
+
+        triSurfaceMesh searchSurf
+        (
+            IOobject
+            (
+                sFeatFileName + ".closeness",
+                runTime.constant(),
+                "triSurface",
+                runTime,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            surf
+        );
+
     // Find close features
 
     // // Dummy trim operation to mark features
@@ -922,72 +886,57 @@ int main(int argc, char *argv[])
     //     )
     // );
 
-    if (args.optionFound("closeness"))
-    {
-        Info<< nl << "Extracting internal and external closeness of surface."
-            << endl;
-
-        // Internal and external closeness
-
-        // Prepare start and end points for intersection tests
-
-        const vectorField& normals = searchSurf.faceNormals();
-
-        scalar span = searchSurf.bounds().mag();
-
-        //args.optionReadIfPresent("closeness", span);
-
-        scalar externalAngleTolerance = 10;
-        scalar externalToleranceCosAngle = Foam::cos
-        (
-            degToRad(180 - externalAngleTolerance)
-        );
-
-        scalar internalAngleTolerance = 45;
-        scalar internalToleranceCosAngle = Foam::cos
-        (
-            degToRad(180 - internalAngleTolerance)
-        );
-
-        Info<< "externalToleranceCosAngle: " << externalToleranceCosAngle << nl
-            << "internalToleranceCosAngle: " << internalToleranceCosAngle
-            << endl;
-
-        // Info<< "span " << span << endl;
-
-        pointField start = searchSurf.faceCentres() - span*normals;
-        pointField end = searchSurf.faceCentres() + span*normals;
-        const pointField& faceCentres = searchSurf.faceCentres();
-
-        List<List<pointIndexHit> > allHitInfo;
-
-        // Find all intersections (in order)
-        searchSurf.findLineAll(start, end, allHitInfo);
-
-        scalarField internalCloseness(start.size(), GREAT);
-        scalarField externalCloseness(start.size(), GREAT);
-
-        forAll(allHitInfo, fI)
+        if (closeness)
         {
-            const List<pointIndexHit>& hitInfo = allHitInfo[fI];
+            Info<< nl << "Extracting internal and external closeness of "
+                << "surface." << endl;
 
-            if (hitInfo.size() < 1)
-            {
-                drawHitProblem(fI, surf, start, faceCentres, end, hitInfo);
+            // Internal and external closeness
 
-                // FatalErrorIn(args.executable())
-                //     << "findLineAll did not hit its own face."
-                //     << exit(FatalError);
-            }
-            else if (hitInfo.size() == 1)
+            // Prepare start and end points for intersection tests
+
+            const vectorField& normals = searchSurf.faceNormals();
+
+            scalar span = searchSurf.bounds().mag();
+
+            scalar externalAngleTolerance = 10;
+            scalar externalToleranceCosAngle =
+                Foam::cos
+                (
+                    degToRad(180 - externalAngleTolerance)
+                );
+
+            scalar internalAngleTolerance = 45;
+            scalar internalToleranceCosAngle =
+                Foam::cos
+                (
+                    degToRad(180 - internalAngleTolerance)
+                );
+
+            Info<< "externalToleranceCosAngle: " << externalToleranceCosAngle
+                << nl
+                << "internalToleranceCosAngle: " << internalToleranceCosAngle
+                << endl;
+
+            // Info<< "span " << span << endl;
+
+            pointField start = searchSurf.faceCentres() - span*normals;
+            pointField end = searchSurf.faceCentres() + span*normals;
+            const pointField& faceCentres = searchSurf.faceCentres();
+
+            List<List<pointIndexHit> > allHitInfo;
+
+            // Find all intersections (in order)
+            searchSurf.findLineAll(start, end, allHitInfo);
+
+            scalarField internalCloseness(start.size(), GREAT);
+            scalarField externalCloseness(start.size(), GREAT);
+
+            forAll(allHitInfo, fI)
             {
-                if (!hitInfo[0].hit())
-                {
-                    // FatalErrorIn(args.executable())
-                    //     << "findLineAll did not hit any face."
-                    //     << exit(FatalError);
-                }
-                else if (hitInfo[0].index() != fI)
+                const List<pointIndexHit>& hitInfo = allHitInfo[fI];
+
+                if (hitInfo.size() < 1)
                 {
                     drawHitProblem(fI, surf, start, faceCentres, end, hitInfo);
 
@@ -995,283 +944,339 @@ int main(int argc, char *argv[])
                     //     << "findLineAll did not hit its own face."
                     //     << exit(FatalError);
                 }
-            }
-            else
-            {
-                label ownHitI = -1;
-
-                forAll(hitInfo, hI)
+                else if (hitInfo.size() == 1)
                 {
-                    // Find the hit on the triangle that launched the ray
-
-                    if (hitInfo[hI].index() == fI)
+                    if (!hitInfo[0].hit())
                     {
-                        ownHitI = hI;
-
-                        break;
+                        // FatalErrorIn(args.executable())
+                        //     << "findLineAll did not hit any face."
+                        //     << exit(FatalError);
                     }
-                }
-
-                if (ownHitI < 0)
-                {
-                    drawHitProblem(fI, surf, start, faceCentres, end, hitInfo);
-
-                    // FatalErrorIn(args.executable())
-                    //     << "findLineAll did not hit its own face."
-                    //     << exit(FatalError);
-                }
-                else if (ownHitI == 0)
-                {
-                    // There are no internal hits, the first hit is the closest
-                    // external hit
-
-                    if
-                    (
-                        (normals[fI] & normals[hitInfo[ownHitI + 1].index()])
-                      < externalToleranceCosAngle
-                    )
+                    else if (hitInfo[0].index() != fI)
                     {
-                        externalCloseness[fI] = mag
+                        drawHitProblem
                         (
-                            faceCentres[fI] - hitInfo[ownHitI + 1].hitPoint()
+                            fI,
+                            surf,
+                            start,
+                            faceCentres,
+                            end,
+                            hitInfo
                         );
-                    }
-                }
-                else if (ownHitI == hitInfo.size() - 1)
-                {
-                    // There are no external hits, the last but one hit is the
-                    // closest internal hit
 
-                    if
-                    (
-                        (normals[fI] & normals[hitInfo[ownHitI - 1].index()])
-                      < internalToleranceCosAngle
-                    )
-                    {
-                        internalCloseness[fI] = mag
-                        (
-                            faceCentres[fI] - hitInfo[ownHitI - 1].hitPoint()
-                        );
+                        // FatalErrorIn(args.executable())
+                        //     << "findLineAll did not hit its own face."
+                        //     << exit(FatalError);
                     }
                 }
                 else
                 {
-                    if
-                    (
-                        (normals[fI] & normals[hitInfo[ownHitI + 1].index()])
-                      < externalToleranceCosAngle
-                    )
+                    label ownHitI = -1;
+
+                    forAll(hitInfo, hI)
                     {
-                        externalCloseness[fI] = mag
-                        (
-                            faceCentres[fI] - hitInfo[ownHitI + 1].hitPoint()
-                        );
+                        // Find the hit on the triangle that launched the ray
+
+                        if (hitInfo[hI].index() == fI)
+                        {
+                            ownHitI = hI;
+
+                            break;
+                        }
                     }
 
-                    if
-                    (
-                        (normals[fI] & normals[hitInfo[ownHitI - 1].index()])
-                      < internalToleranceCosAngle
-                    )
+                    if (ownHitI < 0)
                     {
-                        internalCloseness[fI] = mag
+                        drawHitProblem
                         (
-                            faceCentres[fI] - hitInfo[ownHitI - 1].hitPoint()
+                            fI,
+                            surf,
+                            start,
+                            faceCentres,
+                            end,
+                            hitInfo
                         );
+
+                        // FatalErrorIn(args.executable())
+                        //     << "findLineAll did not hit its own face."
+                        //     << exit(FatalError);
+                    }
+                    else if (ownHitI == 0)
+                    {
+                        // There are no internal hits, the first hit is the
+                        // closest external hit
+
+                        if
+                        (
+                            (
+                                normals[fI]
+                              & normals[hitInfo[ownHitI + 1].index()]
+                            )
+                          < externalToleranceCosAngle
+                        )
+                        {
+                            externalCloseness[fI] =
+                                mag
+                                (
+                                    faceCentres[fI]
+                                  - hitInfo[ownHitI + 1].hitPoint()
+                                );
+                        }
+                    }
+                    else if (ownHitI == hitInfo.size() - 1)
+                    {
+                        // There are no external hits, the last but one hit is
+                        // the closest internal hit
+
+                        if
+                        (
+                            (
+                                normals[fI]
+                              & normals[hitInfo[ownHitI - 1].index()]
+                            )
+                          < internalToleranceCosAngle
+                        )
+                        {
+                            internalCloseness[fI] =
+                                mag
+                                (
+                                    faceCentres[fI]
+                                  - hitInfo[ownHitI - 1].hitPoint()
+                                );
+                        }
+                    }
+                    else
+                    {
+                        if
+                        (
+                            (
+                                normals[fI]
+                              & normals[hitInfo[ownHitI + 1].index()]
+                            )
+                          < externalToleranceCosAngle
+                        )
+                        {
+                            externalCloseness[fI] =
+                                mag
+                                (
+                                    faceCentres[fI]
+                                  - hitInfo[ownHitI + 1].hitPoint()
+                                );
+                        }
+
+                        if
+                        (
+                            (
+                                normals[fI]
+                              & normals[hitInfo[ownHitI - 1].index()]
+                            )
+                          < internalToleranceCosAngle
+                        )
+                        {
+                            internalCloseness[fI] =
+                                mag
+                                (
+                                    faceCentres[fI]
+                                  - hitInfo[ownHitI - 1].hitPoint()
+                                );
+                        }
                     }
                 }
             }
-        }
 
-        triSurfaceScalarField internalClosenessField
-        (
-            IOobject
+            triSurfaceScalarField internalClosenessField
             (
-                sFeatFileName + ".internalCloseness",
-                runTime.constant(),
-                "triSurface",
-                runTime,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            surf,
-            dimLength,
-            internalCloseness
-        );
-
-        internalClosenessField.write();
-
-        triSurfaceScalarField externalClosenessField
-        (
-            IOobject
-            (
-                sFeatFileName + ".externalCloseness",
-                runTime.constant(),
-                "triSurface",
-                runTime,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            surf,
-            dimLength,
-            externalCloseness
-        );
-
-        externalClosenessField.write();
-
-        if (writeVTK)
-        {
-            vtkSurfaceWriter().write
-            (
-                runTime.constant()/"triSurface",    // outputDir
-                sFeatFileName,                      // surfaceName
-                surf.points(),
-                faces,
-                "internalCloseness",                // fieldName
-                internalCloseness,
-                false,                              // isNodeValues
-                true                                // verbose
+                IOobject
+                (
+                    sFeatFileName + ".internalCloseness",
+                    runTime.constant(),
+                    "triSurface",
+                    runTime,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                surf,
+                dimLength,
+                internalCloseness
             );
 
-            vtkSurfaceWriter().write
+            internalClosenessField.write();
+
+            triSurfaceScalarField externalClosenessField
             (
-                runTime.constant()/"triSurface",    // outputDir
-                sFeatFileName,                      // surfaceName
-                surf.points(),
-                faces,
-                "externalCloseness",                // fieldName
-                externalCloseness,
-                false,                              // isNodeValues
-                true                                // verbose
+                IOobject
+                (
+                    sFeatFileName + ".externalCloseness",
+                    runTime.constant(),
+                    "triSurface",
+                    runTime,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                surf,
+                dimLength,
+                externalCloseness
             );
+
+            externalClosenessField.write();
+
+            if (writeVTK)
+            {
+                vtkSurfaceWriter().write
+                (
+                    runTime.constant()/"triSurface",    // outputDir
+                    sFeatFileName,                      // surfaceName
+                    surf.points(),
+                    faces,
+                    "internalCloseness",                // fieldName
+                    internalCloseness,
+                    false,                              // isNodeValues
+                    true                                // verbose
+                );
+
+                vtkSurfaceWriter().write
+                (
+                    runTime.constant()/"triSurface",    // outputDir
+                    sFeatFileName,                      // surfaceName
+                    surf.points(),
+                    faces,
+                    "externalCloseness",                // fieldName
+                    externalCloseness,
+                    false,                              // isNodeValues
+                    true                                // verbose
+                );
+            }
         }
-    }
 
 
 #ifdef ENABLE_CURVATURE
-    if (args.optionFound("curvature"))
-    {
-        Info<< nl << "Extracting curvature of surface at the points." << endl;
+        if (curvature)
+        {
+            Info<< nl << "Extracting curvature of surface at the points."
+                << endl;
 
-        scalarField k = calcCurvature(surf);
+            scalarField k = calcCurvature(surf);
 
         // Modify the curvature values on feature edges and points to be zero.
 
-    //    forAll(newSet.featureEdges(), fEI)
-    //    {
-    //        const edge& e = surf.edges()[newSet.featureEdges()[fEI]];
-    //
-    //        k[surf.meshPoints()[e.start()]] = 0.0;
-    //        k[surf.meshPoints()[e.end()]] = 0.0;
-    //    }
+        //    forAll(newSet.featureEdges(), fEI)
+        //    {
+        //        const edge& e = surf.edges()[newSet.featureEdges()[fEI]];
+        //
+        //        k[surf.meshPoints()[e.start()]] = 0.0;
+        //        k[surf.meshPoints()[e.end()]] = 0.0;
+        //    }
 
-        triSurfacePointScalarField kField
-        (
-            IOobject
+            triSurfacePointScalarField kField
             (
-                sFeatFileName + ".curvature",
-                runTime.constant(),
-                "triSurface",
-                runTime,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            surf,
-            dimLength,
-            k
-        );
-
-        kField.write();
-
-        if (writeVTK)
-        {
-            vtkSurfaceWriter().write
-            (
-                runTime.constant()/"triSurface",    // outputDir
-                sFeatFileName,                      // surfaceName
-                surf.points(),
-                faces,
-                "curvature",                        // fieldName
-                k,
-                true,                               // isNodeValues
-                true                                // verbose
+                IOobject
+                (
+                    sFeatFileName + ".curvature",
+                    runTime.constant(),
+                    "triSurface",
+                    runTime,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                surf,
+                dimLength,
+                k
             );
+
+            kField.write();
+
+            if (writeVTK)
+            {
+                vtkSurfaceWriter().write
+                (
+                    runTime.constant()/"triSurface",    // outputDir
+                    sFeatFileName,                      // surfaceName
+                    surf.points(),
+                    faces,
+                    "curvature",                        // fieldName
+                    k,
+                    true,                               // isNodeValues
+                    true                                // verbose
+                );
+            }
         }
-    }
 #endif
 
 
-    if (args.optionFound("featureProximity"))
-    {
-        Info<< nl << "Extracting proximity of close feature points and edges "
-            << "to the surface" << endl;
-
-        const scalar searchDistance =
-            args.optionRead<scalar>("featureProximity");
-
-        const scalar radiusSqr = sqr(searchDistance);
-
-        scalarField featureProximity(surf.size(), searchDistance);
-
-        forAll(surf, fI)
+        if (featureProximity)
         {
-            const triPointRef& tri = surf[fI].tri(surf.points());
-            const point& triCentre = tri.circumCentre();
+            Info<< nl << "Extracting proximity of close feature points and "
+                << "edges to the surface" << endl;
 
-            List<pointIndexHit> hitList;
+            const scalar searchDistance =
+                readScalar(surfaceDict.lookup("maxFeatureProximity"));
 
-            feMesh.allNearestFeatureEdges(triCentre, radiusSqr, hitList);
+            const scalar radiusSqr = sqr(searchDistance);
 
-            featureProximity[fI] =
-                calcProximityOfFeatureEdges
-                (
-                    feMesh,
-                    hitList,
-                    featureProximity[fI]
-                );
+            scalarField featureProximity(surf.size(), searchDistance);
 
-            feMesh.allNearestFeaturePoints(triCentre, radiusSqr, hitList);
+            forAll(surf, fI)
+            {
+                const triPointRef& tri = surf[fI].tri(surf.points());
+                const point& triCentre = tri.circumCentre();
 
-            featureProximity[fI] =
-                calcProximityOfFeaturePoints
-                (
-                    hitList,
-                    featureProximity[fI]
-                );
-        }
+                List<pointIndexHit> hitList;
 
-        triSurfaceScalarField featureProximityField
-        (
-            IOobject
+                feMesh.allNearestFeatureEdges(triCentre, radiusSqr, hitList);
+
+                featureProximity[fI] =
+                    calcProximityOfFeatureEdges
+                    (
+                        feMesh,
+                        hitList,
+                        featureProximity[fI]
+                    );
+
+                feMesh.allNearestFeaturePoints(triCentre, radiusSqr, hitList);
+
+                featureProximity[fI] =
+                    calcProximityOfFeaturePoints
+                    (
+                        hitList,
+                        featureProximity[fI]
+                    );
+            }
+
+            triSurfaceScalarField featureProximityField
             (
-                sFeatFileName + ".featureProximity",
-                runTime.constant(),
-                "triSurface",
-                runTime,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            surf,
-            dimLength,
-            featureProximity
-        );
-
-        featureProximityField.write();
-
-        if (writeVTK)
-        {
-            vtkSurfaceWriter().write
-            (
-                runTime.constant()/"triSurface",    // outputDir
-                sFeatFileName,                      // surfaceName
-                surf.points(),
-                faces,
-                "featureProximity",                // fieldName
-                featureProximity,
-                false,                              // isNodeValues
-                true                                // verbose
+                IOobject
+                (
+                    sFeatFileName + ".featureProximity",
+                    runTime.constant(),
+                    "triSurface",
+                    runTime,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                surf,
+                dimLength,
+                featureProximity
             );
+
+            featureProximityField.write();
+
+            if (writeVTK)
+            {
+                vtkSurfaceWriter().write
+                (
+                    runTime.constant()/"triSurface",    // outputDir
+                    sFeatFileName,                      // surfaceName
+                    surf.points(),
+                    faces,
+                    "featureProximity",                 // fieldName
+                    featureProximity,
+                    false,                              // isNodeValues
+                    true                                // verbose
+                );
+            }
         }
+
+        Info<< endl;
     }
 
     Info<< "End\n" << endl;
diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractDict b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractDict
new file mode 100644
index 0000000000..ef97d13889
--- /dev/null
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractDict
@@ -0,0 +1,86 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      surfaceFeatureExtractDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+surface1.stl
+{
+    // extractFromFile || extractFromSurface
+    extractionMethod    extractFromFile;
+
+    extractFromFile
+    {
+        // Load from an existing feature edge file
+        featureEdgeFile "constant/triSurface/featureEdges.nas";
+    }
+
+    trimFeatures
+    {
+        // Remove features with fewer than the specified number of edges
+        minElem         0;
+
+        // Remove features shorter than the specified cumulative length
+        minLen          0.0;
+    }
+
+    subsetFeatures
+    {
+        // Use a plane to select feature edges
+        // (normal)(basePoint)
+        plane               (1 0 0)(0 0 0);
+
+        // Select feature edges using a box
+        // (minPt)(maxPt)
+        insideBox           (0 0 0)(1 1 1);
+        outsideBox          (0 0 0)(1 1 1);
+
+        // Remove any non-manifold (open or > 2 connected faces) edges
+        manifoldEdges       no;
+    }
+
+    // Output the curvature of the surface
+    curvature               no;
+
+    // Output the proximity of feature points and edges to each other
+    featureProximity        no;
+    // The maximum search distance to use when looking for other feature
+    // points and edges
+    maxFeatureProximity     1;
+
+    // Out put the closeness of surface elements to other surface elements.
+    closeness               no;
+
+    // Write options
+    writeVTK                no;
+    writeObj                yes;
+    writeFeatureEdgeMesh    no;
+}
+
+
+surface2.nas
+{
+    extractionMethod    extractFromSurface;
+
+    extractFromSurface
+    {
+        // Mark edges whose adjacent surface normals are at an angle less
+        // than includedAngle as features
+        // - 0  : selects no edges
+        // - 180: selects all edges
+        includedAngle   120;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C
index 024f7d27ec..f16dd2595c 100644
--- a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C
+++ b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C
@@ -2832,7 +2832,7 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findNearest
 ) const
 {
     label nearestShapeI = -1;
-    point nearestPoint;
+    point nearestPoint = vector::zero;
 
     if (nodes_.size())
     {
@@ -2847,10 +2847,6 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findNearest
             nearestPoint
         );
     }
-    else
-    {
-        nearestPoint = vector::zero;
-    }
 
     return pointIndexHit(nearestShapeI != -1, nearestPoint, nearestShapeI);
 }
diff --git a/src/OpenFOAM/fields/UniformDimensionedFields/UniformDimensionedField.C b/src/OpenFOAM/fields/UniformDimensionedFields/UniformDimensionedField.C
index 54db26a30c..2879e28f10 100644
--- a/src/OpenFOAM/fields/UniformDimensionedFields/UniformDimensionedField.C
+++ b/src/OpenFOAM/fields/UniformDimensionedFields/UniformDimensionedField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,7 +36,22 @@ Foam::UniformDimensionedField<Type>::UniformDimensionedField
 :
     regIOobject(io),
     dimensioned<Type>(dt)
-{}
+{
+    // Read value
+    if
+    (
+        (
+            io.readOpt() == IOobject::MUST_READ
+         || io.readOpt() == IOobject::MUST_READ_IF_MODIFIED
+        )
+     || (io.readOpt() == IOobject::READ_IF_PRESENT && headerOk())
+    )
+    {
+        dictionary dict(readStream(typeName));
+        this->dimensions().reset(dict.lookup("dimensions"));
+        dict.lookup("value") >> this->value();
+    }
+}
 
 
 template<class Type>
diff --git a/src/OpenFOAM/fields/UniformDimensionedFields/UniformDimensionedField.H b/src/OpenFOAM/fields/UniformDimensionedFields/UniformDimensionedField.H
index 62461df280..bab30a55ff 100644
--- a/src/OpenFOAM/fields/UniformDimensionedFields/UniformDimensionedField.H
+++ b/src/OpenFOAM/fields/UniformDimensionedFields/UniformDimensionedField.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -64,7 +64,7 @@ public:
 
     // Constructors
 
-        //- Construct from components
+        //- Construct from components. Either reads or uses supplied value.
         UniformDimensionedField(const IOobject&, const dimensioned<Type>&);
 
         //- Construct as copy
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C
index 560010c31a..b577a17177 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -1745,6 +1745,7 @@ void Foam::hexRef8::setInstance(const fileName& inst)
 
     cellLevel_.instance() = inst;
     pointLevel_.instance() = inst;
+    level0Edge_.instance() = inst;
     history_.instance() = inst;
 }
 
@@ -1752,7 +1753,7 @@ void Foam::hexRef8::setInstance(const fileName& inst)
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from mesh, read refinement data
-Foam::hexRef8::hexRef8(const polyMesh& mesh)
+Foam::hexRef8::hexRef8(const polyMesh& mesh, const bool readHistory)
 :
     mesh_(mesh),
     cellLevel_
@@ -1781,7 +1782,19 @@ Foam::hexRef8::hexRef8(const polyMesh& mesh)
         ),
         labelList(mesh_.nPoints(), 0)
     ),
-    level0Edge_(getLevel0EdgeLength()),
+    level0Edge_
+    (
+        IOobject
+        (
+            "level0Edge",
+            mesh_.facesInstance(),
+            polyMesh::meshSubDir,
+            mesh_,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        dimensionedScalar("level0Edge", dimLength, getLevel0EdgeLength())
+    ),
     history_
     (
         IOobject
@@ -1790,15 +1803,24 @@ Foam::hexRef8::hexRef8(const polyMesh& mesh)
             mesh_.facesInstance(),
             polyMesh::meshSubDir,
             mesh_,
-            IOobject::READ_IF_PRESENT,
+            IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh_.nCells()    // All cells visible if could not be read
+        (readHistory ? mesh_.nCells() : 0)  // All cells visible if not be read
     ),
     faceRemover_(mesh_, GREAT),     // merge boundary faces wherever possible
     savedPointLevel_(0),
     savedCellLevel_(0)
 {
+    if (readHistory)
+    {
+        history_.readOpt() = IOobject::READ_IF_PRESENT;
+        if (history_.headerOk())
+        {
+            history_.read();
+        }
+    }
+
     if (history_.active() && history_.visibleCells().size() != mesh_.nCells())
     {
         FatalErrorIn
@@ -1849,7 +1871,8 @@ Foam::hexRef8::hexRef8
     const polyMesh& mesh,
     const labelList& cellLevel,
     const labelList& pointLevel,
-    const refinementHistory& history
+    const refinementHistory& history,
+    const scalar level0Edge
 )
 :
     mesh_(mesh),
@@ -1879,7 +1902,24 @@ Foam::hexRef8::hexRef8
         ),
         pointLevel
     ),
-    level0Edge_(getLevel0EdgeLength()),
+    level0Edge_
+    (
+        IOobject
+        (
+            "level0Edge",
+            mesh_.facesInstance(),
+            polyMesh::meshSubDir,
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        dimensionedScalar
+        (
+            "level0Edge",
+            dimLength,
+            (level0Edge >= 0 ? level0Edge : getLevel0EdgeLength())
+        )
+    ),
     history_
     (
         IOobject
@@ -1945,7 +1985,8 @@ Foam::hexRef8::hexRef8
 (
     const polyMesh& mesh,
     const labelList& cellLevel,
-    const labelList& pointLevel
+    const labelList& pointLevel,
+    const scalar level0Edge
 )
 :
     mesh_(mesh),
@@ -1975,7 +2016,24 @@ Foam::hexRef8::hexRef8
         ),
         pointLevel
     ),
-    level0Edge_(getLevel0EdgeLength()),
+    level0Edge_
+    (
+        IOobject
+        (
+            "level0Edge",
+            mesh_.facesInstance(),
+            polyMesh::meshSubDir,
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        dimensionedScalar
+        (
+            "level0Edge",
+            dimLength,
+            (level0Edge >= 0 ? level0Edge : getLevel0EdgeLength())
+        )
+    ),
     history_
     (
         IOobject
@@ -2621,7 +2679,7 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
             << "Value should be >= 1" << exit(FatalError);
     }
 
-    const scalar level0Size = 2*maxFaceDiff*level0Edge_;
+    const scalar level0Size = 2*maxFaceDiff*level0EdgeLength();
 
 
     // Bit tricky. Say we want a distance of three cells between two
@@ -5514,7 +5572,10 @@ void Foam::hexRef8::setUnrefinement
 // Write refinement to polyMesh directory.
 bool Foam::hexRef8::write() const
 {
-    bool writeOk = cellLevel_.write() && pointLevel_.write();
+    bool writeOk =
+        cellLevel_.write()
+     && pointLevel_.write()
+     && level0Edge_.write();
 
     if (history_.active())
     {
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.H
index 36980f7320..503412dd38 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.H
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,6 +43,7 @@ SourceFiles
 #include "removeFaces.H"
 #include "refinementHistory.H"
 #include "PackedBoolList.H"
+#include "uniformDimensionedFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -74,7 +75,7 @@ class hexRef8
         labelIOList pointLevel_;
 
         //- Typical edge length between unrefined points
-        const scalar level0Edge_;
+        uniformDimensionedScalarField level0Edge_;
 
         //- Refinement history
         refinementHistory history_;
@@ -319,24 +320,29 @@ public:
     // Constructors
 
         //- Construct from mesh, read_if_present refinement data
-        //  (from write below)
-        hexRef8(const polyMesh& mesh);
+        //  (from write below). If readHistory is true does read_if_present
+        //  of refinement history. If false clears all history
+        hexRef8(const polyMesh& mesh, const bool readHistory = true);
 
-        //- Construct from mesh and un/refinement data.
+        //- Construct from mesh and un/refinement data and optional size of
+        //  starting cells
         hexRef8
         (
             const polyMesh& mesh,
             const labelList& cellLevel,
             const labelList& pointLevel,
-            const refinementHistory& history
+            const refinementHistory& history,
+            const scalar level0Edge = -GREAT
         );
 
-        //- Construct from mesh and refinement data.
+        //- Construct from mesh and refinement data and optional size of
+        //  starting cells
         hexRef8
         (
             const polyMesh& mesh,
             const labelList& cellLevel,
-            const labelList& pointLevel
+            const labelList& pointLevel,
+            const scalar level0Edge = -GREAT
         );
 
 
@@ -367,7 +373,7 @@ public:
             //- Typical edge length between unrefined points
             scalar level0EdgeLength() const
             {
-                return level0Edge_;
+                return level0Edge_.value();
             }
 
         // Refinement
diff --git a/src/edgeMesh/edgeFormats/nas/NASedgeFormat.C b/src/edgeMesh/edgeFormats/nas/NASedgeFormat.C
index 93196c0e7b..c043ae3c0b 100644
--- a/src/edgeMesh/edgeFormats/nas/NASedgeFormat.C
+++ b/src/edgeMesh/edgeFormats/nas/NASedgeFormat.C
@@ -113,6 +113,17 @@ bool Foam::fileFormats::NASedgeFormat::read
             // discard groupID
             dynEdges.append(e);
         }
+        else if (cmd == "PLOTEL")
+        {
+            edge e;
+
+            // label groupId = readLabel(IStringStream(line.substr(16,8))());
+            e[0] = readLabel(IStringStream(line.substr(16,8))());
+            e[1] = readLabel(IStringStream(line.substr(24,8))());
+
+            // discard groupID
+            dynEdges.append(e);
+        }
         else if (cmd == "GRID")
         {
             label index = readLabel(IStringStream(line.substr(8,8))());
diff --git a/src/edgeMesh/edgeMesh.C b/src/edgeMesh/edgeMesh.C
index 113112fa28..bc63e3d117 100644
--- a/src/edgeMesh/edgeMesh.C
+++ b/src/edgeMesh/edgeMesh.C
@@ -28,6 +28,7 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "addToMemberFunctionSelectionTable.H"
 #include "ListOps.H"
+#include "EdgeMap.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -327,10 +328,7 @@ void Foam::edgeMesh::mergePoints(const scalar mergeDist)
         }
 
         // Compact using a hashtable and commutative hash of edge.
-        HashTable<label, edge, Hash<edge> > edgeToLabel
-        (
-            2*edges_.size()
-        );
+        EdgeMap<label> edgeToLabel(2*edges_.size());
 
         label newEdgeI = 0;
 
@@ -349,13 +347,7 @@ void Foam::edgeMesh::mergePoints(const scalar mergeDist)
 
         edges_.setSize(newEdgeI);
 
-        for
-        (
-            HashTable<label, edge, Hash<edge> >::const_iterator iter =
-                edgeToLabel.begin();
-            iter != edgeToLabel.end();
-            ++iter
-        )
+        forAllConstIter(EdgeMap<label>, edgeToLabel, iter)
         {
             edges_[iter()] = iter.key();
         }
@@ -363,4 +355,35 @@ void Foam::edgeMesh::mergePoints(const scalar mergeDist)
 }
 
 
+void Foam::edgeMesh::mergeEdges()
+{
+    EdgeMap<label> existingEdges(2*edges_.size());
+
+    label curEdgeI = 0;
+    forAll(edges_, edgeI)
+    {
+        const edge& e = edges_[edgeI];
+
+        if (existingEdges.insert(e, curEdgeI))
+        {
+            curEdgeI++;
+        }
+    }
+
+    if (debug)
+    {
+        Info<< "Merging duplicate edges: "
+            << edges_.size() - existingEdges.size()
+            << " edges will be deleted." << endl;
+    }
+
+    edges_.setSize(existingEdges.size());
+
+    forAllConstIter(EdgeMap<label>, existingEdges, iter)
+    {
+        edges_[iter()] = iter.key();
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/src/edgeMesh/edgeMesh.H b/src/edgeMesh/edgeMesh.H
index baa6142815..66574e9fac 100644
--- a/src/edgeMesh/edgeMesh.H
+++ b/src/edgeMesh/edgeMesh.H
@@ -246,6 +246,9 @@ public:
         //- Merge common points (points within mergeDist)
         void mergePoints(const scalar mergeDist);
 
+        //- Merge similar edges
+        void mergeEdges();
+
 
     // Write
 
diff --git a/src/fieldSources/basicSource/explicitSource/ExplicitSource.H b/src/fieldSources/basicSource/explicitSource/ExplicitSource.H
index e53e6004be..95e752d17b 100644
--- a/src/fieldSources/basicSource/explicitSource/ExplicitSource.H
+++ b/src/fieldSources/basicSource/explicitSource/ExplicitSource.H
@@ -39,7 +39,7 @@ Description
             }
         }
 
-    If volumeMode = 
+    If volumeMode =
         - absolute: values are given as <quantity>
         - specific: values are given as <quantity>/m3
 
diff --git a/src/fieldSources/basicSource/pressureGradientExplicitSource/pressureGradientExplicitSource.H b/src/fieldSources/basicSource/pressureGradientExplicitSource/pressureGradientExplicitSource.H
index bf760254af..48b4bf2c24 100644
--- a/src/fieldSources/basicSource/pressureGradientExplicitSource/pressureGradientExplicitSource.H
+++ b/src/fieldSources/basicSource/pressureGradientExplicitSource/pressureGradientExplicitSource.H
@@ -38,7 +38,7 @@ Description
             gradPini    gradPini [0 2 -2 0 0] 0; // initial pressure gradient
             flowDir     (1 0 0);    // flow direction
         }
-        
+
 
 SourceFiles
     pressureGradientExplicitSource.C
diff --git a/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.H b/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.H
index ae2d8a2a47..946c752999 100644
--- a/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.H
+++ b/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.H
@@ -202,7 +202,7 @@ protected:
         //- Maximum radius
         scalar rMax_;
 
-        
+
     // Protected Member Functions
 
         //- Check data
@@ -247,7 +247,7 @@ public:
 
     // Constructors
 
-        
+
         //- Construct from components
         rotorDiskSource
         (
diff --git a/src/fileFormats/vtk/vtkUnstructuredReader.C b/src/fileFormats/vtk/vtkUnstructuredReader.C
index 3d80e1ccec..049f68f89b 100644
--- a/src/fileFormats/vtk/vtkUnstructuredReader.C
+++ b/src/fileFormats/vtk/vtkUnstructuredReader.C
@@ -704,6 +704,7 @@ void Foam::vtkUnstructuredReader::read(ISstream& inFile)
             labelList faceVerts;
             readBlock(inFile, nNumbers, faceVerts);
             faces_.setSize(nFaces);
+            faceMap_.setSize(nFaces);
             label elemI = 0;
             forAll(faces_, faceI)
             {
diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index 3652bea89a..0e2de76c07 100644
--- a/src/finiteVolume/Make/files
+++ b/src/finiteVolume/Make/files
@@ -124,6 +124,7 @@ $(derivedFvPatchFields)/mappedFixedPushedInternalValue/mappedFixedPushedInternal
 $(derivedFvPatchFields)/mappedFixedValue/mappedFixedValueFvPatchFields.C
 $(derivedFvPatchFields)/mappedVelocityFluxFixedValue/mappedVelocityFluxFixedValueFvPatchField.C
 $(derivedFvPatchFields)/mappedFlowRate/mappedFlowRateFvPatchVectorField.C
+$(derivedFvPatchFields)/fixedMean/fixedMeanFvPatchFields.C
 $(derivedFvPatchFields)/fan/fanFvPatchFields.C
 $(derivedFvPatchFields)/fanPressure/fanPressureFvPatchScalarField.C
 $(derivedFvPatchFields)/buoyantPressure/buoyantPressureFvPatchScalarField.C
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchField.C
new file mode 100644
index 0000000000..e59f98910e
--- /dev/null
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchField.C
@@ -0,0 +1,142 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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 "fixedMeanFvPatchField.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+fixedMeanFvPatchField<Type>::fixedMeanFvPatchField
+(
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF
+)
+:
+    fixedValueFvPatchField<Type>(p, iF),
+    meanValue_(pTraits<Type>::zero)
+{}
+
+
+template<class Type>
+fixedMeanFvPatchField<Type>::fixedMeanFvPatchField
+(
+    const fixedMeanFvPatchField<Type>& ptf,
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchField<Type>(ptf, p, iF, mapper),
+    meanValue_(ptf.meanValue_)
+{}
+
+
+template<class Type>
+fixedMeanFvPatchField<Type>::fixedMeanFvPatchField
+(
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchField<Type>(p, iF, dict),
+    meanValue_(pTraits<Type>(dict.lookup("meanValue")))
+{}
+
+
+template<class Type>
+fixedMeanFvPatchField<Type>::fixedMeanFvPatchField
+(
+    const fixedMeanFvPatchField<Type>& ptf
+)
+:
+    fixedValueFvPatchField<Type>(ptf),
+    meanValue_(ptf.meanValue_)
+{}
+
+
+template<class Type>
+fixedMeanFvPatchField<Type>::fixedMeanFvPatchField
+(
+    const fixedMeanFvPatchField<Type>& ptf,
+    const DimensionedField<Type, volMesh>& iF
+)
+:
+    fixedValueFvPatchField<Type>(ptf, iF),
+    meanValue_(ptf.meanValue_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void fixedMeanFvPatchField<Type>::updateCoeffs()
+{
+    if (this->updated())
+    {
+        return;
+    }
+
+    Field<Type> newValues(this->patchInternalField());
+
+    Type meanValuePsi =
+        gSum(this->patch().magSf()*newValues)
+       /gSum(this->patch().magSf());
+
+    if (mag(meanValue_) > SMALL && mag(meanValuePsi)/mag(meanValue_) > 0.5)
+    {
+        newValues *= mag(meanValue_)/mag(meanValuePsi);
+    }
+    else
+    {
+        newValues += (meanValue_ - meanValuePsi);
+    }
+
+    this->operator==(newValues);
+
+    fixedValueFvPatchField<Type>::updateCoeffs();
+}
+
+
+template<class Type>
+void fixedMeanFvPatchField<Type>::write(Ostream& os) const
+{
+    fvPatchField<Type>::write(os);
+    os.writeKeyword("meanValue") << meanValue_ << token::END_STATEMENT << nl;
+    this->writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchField.H
new file mode 100644
index 0000000000..b59969a967
--- /dev/null
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchField.H
@@ -0,0 +1,158 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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::fixedMeanFvPatchField
+
+Description
+    Extrapolates field to the patch using the near-cell values and adjusts
+    the distribution to match the specified meanValue.
+
+SourceFiles
+    fixedMeanFvPatchField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fixedMeanFvPatchField_H
+#define fixedMeanFvPatchField_H
+
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                  Class fixedMeanFvPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class fixedMeanFvPatchField
+:
+    public fixedValueFvPatchField<Type>
+{
+
+protected:
+
+    // Protected data
+
+        //- MeanValue value the field is adjusted to maintain
+        Type meanValue_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("fixedMean");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        fixedMeanFvPatchField
+        (
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        fixedMeanFvPatchField
+        (
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given fixedMeanFvPatchField
+        //  onto a new patch
+        fixedMeanFvPatchField
+        (
+            const fixedMeanFvPatchField<Type>&,
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        fixedMeanFvPatchField
+        (
+            const fixedMeanFvPatchField<Type>&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchField<Type> > clone() const
+        {
+            return tmp<fvPatchField<Type> >
+            (
+                new fixedMeanFvPatchField<Type>(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        fixedMeanFvPatchField
+        (
+            const fixedMeanFvPatchField<Type>&,
+            const DimensionedField<Type, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchField<Type> > clone
+        (
+            const DimensionedField<Type, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchField<Type> >
+            (
+                new fixedMeanFvPatchField<Type>(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "fixedMeanFvPatchField.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchFields.C b/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchFields.C
new file mode 100644
index 0000000000..e8e269c8a5
--- /dev/null
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchFields.C
@@ -0,0 +1,43 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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 "fixedMeanFvPatchFields.H"
+#include "volMesh.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+makePatchFields(fixedMean);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchFields.H b/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchFields.H
new file mode 100644
index 0000000000..64996c3a68
--- /dev/null
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchFields.H
@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fixedMeanFvPatchFields_H
+#define fixedMeanFvPatchFields_H
+
+#include "fixedMeanFvPatchField.H"
+#include "fieldTypes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeFieldTypedefs(fixedMean);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/buildCGALPolyhedron.C b/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchFieldsFwd.H
similarity index 50%
rename from applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/buildCGALPolyhedron.C
rename to src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchFieldsFwd.H
index 63e2b904bf..9f04dc965b 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/CGALPolyhedron/buildCGALPolyhedron.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fixedMean/fixedMeanFvPatchFieldsFwd.H
@@ -23,67 +23,28 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "buildCGALPolyhedron.H"
+#ifndef fixedMeanFvPatchFieldsFwd_H
+#define fixedMeanFvPatchFieldsFwd_H
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+#include "fieldTypes.H"
 
-Foam::buildCGALPolyhedron::buildCGALPolyhedron
-(
-    const Foam::triSurface& surf
-)
-:
-    CGAL::Modifier_base<HalfedgeDS>(),
-    surf_(surf)
-{}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::buildCGALPolyhedron::~buildCGALPolyhedron()
-{}
-
-
-// * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
-
-void Foam::buildCGALPolyhedron::operator()
-(
-    HalfedgeDS& hds
-)
+namespace Foam
 {
-    typedef HalfedgeDS::Traits     Traits;
-    typedef Traits::Point_3        Point;
 
-    // Postcondition: `hds' is a valid polyhedral surface.
-    CGAL::Polyhedron_incremental_builder_3<HalfedgeDS> B(hds, false);
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-    B.begin_surface
-    (
-        surf_.points().size(),      // n points
-        surf_.size(),               // n facets
-        2*surf_.edges().size()      // n halfedges
-    );
+template<class Type> class fixedMeanFvPatchField;
 
-    forAll(surf_.points(), pI)
-    {
-        const Foam::point& p = surf_.points()[pI];
+makePatchTypeFieldTypedefs(fixedMean);
 
-        B.add_vertex(Point(p.x(), p.y(), p.z()));
-    }
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-    forAll(surf_, fI)
-    {
-        B.begin_facet();
+} // End namespace Foam
 
-        B.add_vertex_to_facet(surf_[fI][0]);
-        B.add_vertex_to_facet(surf_[fI][1]);
-        B.add_vertex_to_facet(surf_[fI][2]);
-
-        B.end_facet();
-    }
-
-    B.end_surface();
-}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#endif
 
 // ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.C b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.C
index 4d4f4dc3ca..87d03eade7 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.C
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.C
@@ -181,7 +181,7 @@ void Foam::ParticleErosion<CloudType>::postPatch
 
         // particle direction of travel
         const vector& U = p.U();
-        
+
         // quick reject if particle travelling away from the patch
         if ((-nw & U) < 0)
         {
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/VoidFraction/VoidFraction.H b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/VoidFraction/VoidFraction.H
index 0b15a6f620..26f0f65018 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/VoidFraction/VoidFraction.H
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/VoidFraction/VoidFraction.H
@@ -112,7 +112,7 @@ public:
 
             //- Post-move hook
             virtual void postMove
-            (   
+            (
                 const parcelType& p,
                 const label cellI,
                 const scalar dt
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C
index 4910faa51c..d5e459ddc0 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C
@@ -285,7 +285,7 @@ void Foam::PairCollision<CloudType>::wallInteraction()
                             );
 
                             flatSiteData.append(wSD);
- 
+
                             particleHit = true;
                         }
                     }
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C
index ba8d0f4200..ab2140a4ec 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -76,6 +76,9 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
 Foam::labelList Foam::refinementParameters::findCells(const polyMesh& mesh)
  const
 {
+    // Force calculation of tet-diag decomposition (for use in findCell)
+    (void)mesh.tetBasePtIs();
+
     // Global calculation engine
     globalIndex globalCells(mesh.nCells());
 
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
index 79df04781e..c2565fc670 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -877,49 +877,7 @@ Foam::meshRefinement::meshRefinement
     meshCutter_
     (
         mesh,
-        labelIOList
-        (
-            IOobject
-            (
-                "cellLevel",
-                mesh_.facesInstance(),
-                fvMesh::meshSubDir,
-                mesh,
-                IOobject::READ_IF_PRESENT,
-                IOobject::NO_WRITE,
-                false
-            ),
-            labelList(mesh_.nCells(), 0)
-        ),
-        labelIOList
-        (
-            IOobject
-            (
-                "pointLevel",
-                mesh_.facesInstance(),
-                fvMesh::meshSubDir,
-                mesh,
-                IOobject::READ_IF_PRESENT,
-                IOobject::NO_WRITE,
-                false
-            ),
-            labelList(mesh_.nPoints(), 0)
-        ),
-        refinementHistory
-        (
-            IOobject
-            (
-                "refinementHistory",
-                mesh_.facesInstance(),
-                fvMesh::meshSubDir,
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            List<refinementHistory::splitCell8>(0),
-            labelList(0)
-        )   // no unrefinement
+        false   // do not try to read history.
     ),
     surfaceIndex_
     (
@@ -1819,6 +1777,9 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMeshRegions
     const point& keepPoint
 )
 {
+    // Force calculation of face decomposition (used in findCell)
+    (void)mesh_.tetBasePtIs();
+
     // Determine connected regions. regionSplit is the labelList with the
     // region per cell.
     regionSplit cellRegion(mesh_);
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
index c6877bbd8d..5ec9065955 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -1235,6 +1235,9 @@ void Foam::meshRefinement::findCellZoneInsideWalk
     blockedFace.clear();
 
 
+    // Force calculation of face decomposition (used in findCell)
+    (void)mesh_.tetBasePtIs();
+
     // For all locationSurface find the cell
     forAll(locationSurfaces, i)
     {
diff --git a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C
index e39eae5704..465734ff87 100644
--- a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C
+++ b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C
@@ -33,11 +33,15 @@ License
 #include "OFstream.H"
 #include "IFstream.H"
 #include "unitConversion.H"
+#include "EdgeMap.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 defineTypeNameAndDebug(Foam::surfaceFeatures, 0);
 
+const Foam::scalar Foam::surfaceFeatures::parallelTolerance =
+    sin(degToRad(1.0));
+
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -170,7 +174,10 @@ void Foam::surfaceFeatures::setFromStatus(const List<edgeStatus>& edgeStat)
 
 
 //construct feature points where more than 2 feature edges meet
-void Foam::surfaceFeatures::calcFeatPoints(const List<edgeStatus>& edgeStat)
+void Foam::surfaceFeatures::calcFeatPoints
+(
+    const List<edgeStatus>& edgeStat
+)
 {
     DynamicList<label> featurePoints(surf_.nPoints()/1000);
 
@@ -200,6 +207,60 @@ void Foam::surfaceFeatures::calcFeatPoints(const List<edgeStatus>& edgeStat)
 }
 
 
+void Foam::surfaceFeatures::classifyFeatureAngles
+(
+    const labelListList& edgeFaces,
+    List<edgeStatus>& edgeStat,
+    const scalar minCos
+) const
+{
+    const vectorField& faceNormals = surf_.faceNormals();
+    const pointField& points = surf_.points();
+
+    forAll(edgeFaces, edgeI)
+    {
+        const labelList& eFaces = edgeFaces[edgeI];
+
+        if (eFaces.size() != 2)
+        {
+            // Non-manifold. What to do here? Is region edge? external edge?
+            edgeStat[edgeI] = REGION;
+        }
+        else
+        {
+            label face0 = eFaces[0];
+            label face1 = eFaces[1];
+
+            if (surf_[face0].region() != surf_[face1].region())
+            {
+                edgeStat[edgeI] = REGION;
+            }
+            else
+            {
+                if ((faceNormals[face0] & faceNormals[face1]) < minCos)
+                {
+
+                    // Check if convex or concave by looking at angle
+                    // between face centres and normal
+                    vector f0Tof1 =
+                        surf_[face1].centre(points)
+                      - surf_[face0].centre(points);
+
+                    if ((f0Tof1 & faceNormals[face0]) > 0.0)
+                    {
+                        edgeStat[edgeI] = INTERNAL;
+                    }
+                    else
+                    {
+                        edgeStat[edgeI] = EXTERNAL;
+                    }
+                }
+            }
+        }
+    }
+}
+
+
 // Returns next feature edge connected to pointI with correct value.
 Foam::label Foam::surfaceFeatures::nextFeatEdge
 (
@@ -436,6 +497,80 @@ Foam::surfaceFeatures::surfaceFeatures
 }
 
 
+Foam::surfaceFeatures::surfaceFeatures
+(
+    const triSurface& surf,
+    const pointField& points,
+    const edgeList& edges,
+    const scalar mergeTol
+)
+:
+    surf_(surf),
+    featurePoints_(0),
+    featureEdges_(0),
+    externalStart_(0),
+    internalStart_(0)
+{
+    // Match edge mesh edges with the triSurface edges
+
+    const labelListList& surfEdgeFaces = surf_.edgeFaces();
+    const edgeList& surfEdges = surf_.edges();
+
+    scalar mergeTolSqr = sqr(mergeTol);
+
+    EdgeMap<label> dynFeatEdges(2*edges.size());
+    DynamicList<labelList> dynFeatureEdgeFaces(edges.size());
+
+    labelList edgeLabel;
+
+    nearestFeatEdge
+    (
+        edges,
+        points,
+        mergeTolSqr,
+        edgeLabel     // label of surface edge or -1
+    );
+
+    label count = 0;
+    forAll(edgeLabel, sEdgeI)
+    {
+        const label sEdge = edgeLabel[sEdgeI];
+
+        if (sEdge == -1)
+        {
+            continue;
+        }
+
+        dynFeatEdges.insert(surfEdges[sEdge], count++);
+        dynFeatureEdgeFaces.append(surfEdgeFaces[sEdge]);
+    }
+
+    // Find whether an edge is external or internal
+    List<edgeStatus> edgeStat(dynFeatEdges.size(), NONE);
+
+    classifyFeatureAngles(dynFeatureEdgeFaces, edgeStat, GREAT);
+
+    // Transfer the edge status to a list encompassing all edges in the surface
+    // so that calcFeatPoints can be used.
+    List<edgeStatus> allEdgeStat(surf_.nEdges(), NONE);
+
+    forAll(allEdgeStat, eI)
+    {
+        EdgeMap<label>::const_iterator iter = dynFeatEdges.find(surfEdges[eI]);
+
+        if (iter != dynFeatEdges.end())
+        {
+            allEdgeStat[eI] = edgeStat[iter()];
+        }
+    }
+
+    edgeStat.clear();
+    dynFeatEdges.clear();
+
+    setFromStatus(allEdgeStat);
+}
+
+
 //- Construct as copy
 Foam::surfaceFeatures::surfaceFeatures(const surfaceFeatures& sf)
 :
@@ -496,54 +631,10 @@ void Foam::surfaceFeatures::findFeatures(const scalar includedAngle)
 {
     scalar minCos = Foam::cos(degToRad(180.0 - includedAngle));
 
-    const labelListList& edgeFaces = surf_.edgeFaces();
-    const vectorField& faceNormals = surf_.faceNormals();
-    const pointField& points = surf_.points();
-
     // Per edge whether is feature edge.
     List<edgeStatus> edgeStat(surf_.nEdges(), NONE);
 
-    forAll(edgeFaces, edgeI)
-    {
-        const labelList& eFaces = edgeFaces[edgeI];
-
-        if (eFaces.size() != 2)
-        {
-            // Non-manifold. What to do here? Is region edge? external edge?
-            edgeStat[edgeI] = REGION;
-        }
-        else
-        {
-            label face0 = eFaces[0];
-            label face1 = eFaces[1];
-
-            if (surf_[face0].region() != surf_[face1].region())
-            {
-                edgeStat[edgeI] = REGION;
-            }
-            else
-            {
-                if ((faceNormals[face0] & faceNormals[face1]) < minCos)
-                {
-
-                    // Check if convex or concave by looking at angle
-                    // between face centres and normal
-                    vector f0Tof1 =
-                        surf_[face1].centre(points)
-                        - surf_[face0].centre(points);
-
-                    if ((f0Tof1 & faceNormals[face0]) > 0.0)
-                    {
-                        edgeStat[edgeI] = INTERNAL;
-                    }
-                    else
-                    {
-                        edgeStat[edgeI] = EXTERNAL;
-                    }
-                }
-            }
-        }
-    }
+    classifyFeatureAngles(surf_.edgeFaces(), edgeStat, minCos);
 
     setFromStatus(edgeStat);
 }
@@ -1147,6 +1238,9 @@ void Foam::surfaceFeatures::nearestSurfEdge
 
     const pointField& localPoints = surf_.localPoints();
 
+    treeBoundBox searchDomain(localPoints);
+    searchDomain.inflate(0.1);
+
     indexedOctree<treeDataEdge> ppTree
     (
         treeDataEdge
@@ -1156,7 +1250,7 @@ void Foam::surfaceFeatures::nearestSurfEdge
             localPoints,
             selectedEdges
         ),          // all information needed to do geometric checks
-        treeBoundBox(localPoints),  // overall search domain
+        searchDomain,  // overall search domain
         8,      // maxLevel
         10,     // leafsize
         3.0     // duplicity
@@ -1216,6 +1310,8 @@ void Foam::surfaceFeatures::nearestSurfEdge
     pointOnEdge.setSize(selectedSampleEdges.size());
     pointOnFeature.setSize(selectedSampleEdges.size());
 
+    treeBoundBox searchDomain(surf_.localPoints());
+
     indexedOctree<treeDataEdge> ppTree
     (
         treeDataEdge
@@ -1224,11 +1320,11 @@ void Foam::surfaceFeatures::nearestSurfEdge
             surf_.edges(),
             surf_.localPoints(),
             selectedEdges
-        ),          // all information needed to do geometric checks
-        treeBoundBox(surf_.localPoints()),  // overall search domain
-        8,      // maxLevel
-        10,     // leafsize
-        3.0     // duplicity
+        ),              // all information needed to do geometric checks
+        searchDomain,   // overall search domain
+        8,              // maxLevel
+        10,             // leafsize
+        3.0             // duplicity
     );
 
     forAll(selectedSampleEdges, i)
@@ -1262,6 +1358,67 @@ void Foam::surfaceFeatures::nearestSurfEdge
 }
 
 
+void Foam::surfaceFeatures::nearestFeatEdge
+(
+    const edgeList& edges,
+    const pointField& points,
+    scalar searchSpanSqr,   // Search span
+    labelList& edgeLabel
+) const
+{
+    edgeLabel = labelList(surf_.nEdges(), -1);
+
+    treeBoundBox searchDomain(points);
+    searchDomain.inflate(0.1);
+
+    indexedOctree<treeDataEdge> ppTree
+    (
+        treeDataEdge
+        (
+            false,
+            edges,
+            points,
+            identity(edges.size())
+        ),          // all information needed to do geometric checks
+        searchDomain,  // overall search domain
+        8,      // maxLevel
+        10,     // leafsize
+        3.0     // duplicity
+    );
+
+    const edgeList& surfEdges = surf_.edges();
+    const pointField& surfLocalPoints = surf_.localPoints();
+
+    forAll(surfEdges, edgeI)
+    {
+        const edge& sample = surfEdges[edgeI];
+
+        const point& startPoint = surfLocalPoints[sample.start()];
+        const point& midPoint = sample.centre(surfLocalPoints);
+
+        pointIndexHit infoMid = ppTree.findNearest
+        (
+            midPoint,
+            searchSpanSqr
+        );
+
+        if (infoMid.hit())
+        {
+            const vector surfEdgeDir = midPoint - startPoint;
+
+            const edge& featEdge = edges[infoMid.index()];
+            const vector featEdgeDir = featEdge.vec(points);
+
+            // Check that the edges are nearly parallel
+            if (mag(surfEdgeDir ^ featEdgeDir) < parallelTolerance)
+            {
+                edgeLabel[edgeI] = edgeI;
+            }
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 void Foam::surfaceFeatures::operator=(const surfaceFeatures& rhs)
diff --git a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H
index 7f7101b662..0fca9798a2 100644
--- a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H
+++ b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H
@@ -95,6 +95,11 @@ private:
         {}
     };
 
+    // Static data
+
+        //- Tolerance for determining whether two vectors are parallel
+        static const scalar parallelTolerance;
+
 
     // Private data
 
@@ -130,6 +135,14 @@ private:
         //- Construct feature points where more than 2 feature edges meet
         void calcFeatPoints(const List<edgeStatus>&);
 
+        //- Classify the angles of the feature edges
+        void classifyFeatureAngles
+        (
+            const labelListList& edgeFaces,
+            List<edgeStatus>& edgeStat,
+            const scalar minCos
+        ) const;
+
         //- Choose next unset feature edge.
         label nextFeatEdge
         (
@@ -186,6 +199,15 @@ public:
         //- Construct from file
         surfaceFeatures(const triSurface&, const fileName& fName);
 
+        //- Construct from pointField and edgeList (edgeMesh)
+        surfaceFeatures
+        (
+            const triSurface&,
+            const pointField& points,
+            const edgeList& edges,
+            const scalar mergeTol = 1e-6
+        );
+
         //- Construct as copy
         surfaceFeatures(const surfaceFeatures&);
 
@@ -327,9 +349,8 @@ public:
                 pointField& edgePoint
             ) const;
 
-
             //- Find nearest surface edge (out of selectedEdges) for each
-            // sample edge.
+            //  sample edge.
             //  Sets:
             //  - edgeLabel         : label of surface edge.
             //  - pointOnEdge       : exact position of nearest point on edge.
@@ -347,6 +368,16 @@ public:
                 pointField& pointOnFeature  // point on sample edge
             ) const;
 
+            //- Find nearest feature edge to each surface edge. Uses the
+            //  mid-point of the surface edges.
+            void nearestFeatEdge
+            (
+                const edgeList& edges,
+                const pointField& points,
+                scalar searchSpanSqr,
+                labelList& edgeLabel
+            ) const;
+
 
         // Write
 
diff --git a/src/parallel/decompose/decompositionMethods/structuredDecomp/structuredDecomp.C b/src/parallel/decompose/decompositionMethods/structuredDecomp/structuredDecomp.C
index 072345df2d..1672e53d21 100644
--- a/src/parallel/decompose/decompositionMethods/structuredDecomp/structuredDecomp.C
+++ b/src/parallel/decompose/decompositionMethods/structuredDecomp/structuredDecomp.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -49,11 +49,11 @@ namespace Foam
 Foam::structuredDecomp::structuredDecomp(const dictionary& decompositionDict)
 :
     decompositionMethod(decompositionDict),
-    methodDict_(decompositionDict_.subDict(typeName + "Coeffs"))
+    methodDict_(decompositionDict_.subDict(typeName + "Coeffs")),
+    patches_(methodDict_.lookup("patches"))
 {
     methodDict_.set("numberOfSubdomains", nDomains());
     method_ = decompositionMethod::New(methodDict_);
-    patches_ = wordList(methodDict_.lookup("patches"));
 }
 
 
@@ -72,29 +72,20 @@ Foam::labelList Foam::structuredDecomp::decompose
     const scalarField& cWeights
 )
 {
-    labelList patchIDs(patches_.size());
     const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+    const labelHashSet patchIDs(pbm.patchSet(patches_));
 
     label nFaces = 0;
-    forAll(patches_, i)
+    forAllConstIter(labelHashSet, patchIDs, iter)
     {
-        patchIDs[i] = pbm.findPatchID(patches_[i]);
-
-        if (patchIDs[i] == -1)
-        {
-            FatalErrorIn("structuredDecomp::decompose(..)")
-                << "Cannot find patch " << patches_[i] << endl
-                << "Valid patches are " << pbm.names()
-                << exit(FatalError);
-        }
-        nFaces += pbm[patchIDs[i]].size();
+        nFaces += pbm[iter.key()].size();
     }
 
     // Extract a submesh.
     labelHashSet patchCells(2*nFaces);
-    forAll(patchIDs, i)
+    forAllConstIter(labelHashSet, patchIDs, iter)
     {
-        const labelUList& fc = pbm[patchIDs[i]].faceCells();
+        const labelUList& fc = pbm[iter.key()].faceCells();
         forAll(fc, i)
         {
             patchCells.insert(fc[i]);
@@ -127,9 +118,9 @@ Foam::labelList Foam::structuredDecomp::decompose
     labelList patchFaces(nFaces);
     List<topoDistanceData> patchData(nFaces);
     nFaces = 0;
-    forAll(patchIDs, i)
+    forAllConstIter(labelHashSet, patchIDs, iter)
     {
-        const polyPatch& pp = pbm[patchIDs[i]];
+        const polyPatch& pp = pbm[iter.key()];
         const labelUList& fc = pp.faceCells();
         forAll(fc, i)
         {
diff --git a/src/parallel/decompose/decompositionMethods/structuredDecomp/structuredDecomp.H b/src/parallel/decompose/decompositionMethods/structuredDecomp/structuredDecomp.H
index 259102ac79..a2374f40a0 100644
--- a/src/parallel/decompose/decompositionMethods/structuredDecomp/structuredDecomp.H
+++ b/src/parallel/decompose/decompositionMethods/structuredDecomp/structuredDecomp.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,7 +25,7 @@ Class
     Foam::structuredDecomp
 
 Description
-    Decomposition given using consecutive application of decomposers.
+    Decomposition by walking out decomposition of patch cells mesh.
 
 SourceFiles
     structuredDecomp.C
@@ -52,9 +52,9 @@ class structuredDecomp
 
         dictionary methodDict_;
 
-        autoPtr<decompositionMethod> method_;
+        wordReList patches_;
 
-        wordList patches_;
+        autoPtr<decompositionMethod> method_;
 
 
     // Private Member Functions
diff --git a/src/postProcessing/functionObjects/field/turbulenceFields/postProcessingDict b/src/postProcessing/functionObjects/field/turbulenceFields/postProcessingDict
index 6d4b0409fa..bc35e7dc0f 100644
--- a/src/postProcessing/functionObjects/field/turbulenceFields/postProcessingDict
+++ b/src/postProcessing/functionObjects/field/turbulenceFields/postProcessingDict
@@ -26,7 +26,7 @@ functions
 
         fields
         (
-            R 
+            R
         );
     }
 }
diff --git a/src/postProcessing/functionObjects/utilities/pressureCoefficient/pressureCoefficient.H b/src/postProcessing/functionObjects/utilities/pressureCoefficient/pressureCoefficient.H
index 5f0379aa7f..7d72e3b784 100644
--- a/src/postProcessing/functionObjects/utilities/pressureCoefficient/pressureCoefficient.H
+++ b/src/postProcessing/functionObjects/utilities/pressureCoefficient/pressureCoefficient.H
@@ -28,7 +28,7 @@ Description
     Calculates pressure coefficient, c_p
 
         c_p = p/p_dyn,inf
-    
+
     where:
 
         p_dyn,inf = 0.5*rho*mag(U_inf)^2
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/calculateMeshToMeshAddressing.C b/src/sampling/meshToMeshInterpolation/meshToMesh/calculateMeshToMeshAddressing.C
index 831565bd13..fda5525f09 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/calculateMeshToMeshAddressing.C
+++ b/src/sampling/meshToMeshInterpolation/meshToMesh/calculateMeshToMeshAddressing.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -202,7 +202,7 @@ void Foam::meshToMesh::calcAddressing()
     if (debug)
     {
         Info<< "meshToMesh::calculateAddressing() : "
-            << "finished calculating mesh-to-mesh acell ddressing" << endl;
+            << "finished calculating mesh-to-mesh cell addressing" << endl;
     }
 }
 
diff --git a/src/surfMesh/surfaceFormats/vtk/VTKsurfaceFormat.C b/src/surfMesh/surfaceFormats/vtk/VTKsurfaceFormat.C
index 3dad786857..cc3741da12 100644
--- a/src/surfMesh/surfaceFormats/vtk/VTKsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/vtk/VTKsurfaceFormat.C
@@ -128,6 +128,18 @@ bool Foam::fileFormats::VTKsurfaceFormat<Face>::read
             zones[i] = label(region[i]);
         }
     }
+    else if (reader.cellData().foundObject<scalarIOField>("STLSolidLabeling"))
+    {
+        const scalarIOField& region =
+            reader.cellData().lookupObject<scalarIOField>
+            (
+                "STLSolidLabeling"
+            );
+        forAll(region, i)
+        {
+            zones[i] = label(region[i]);
+        }
+    }
 
     // Create zone names
     const label nZones = max(zones)+1;
diff --git a/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C b/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C
index 49d23cc84a..d1a924ad11 100644
--- a/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C
+++ b/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C
@@ -284,24 +284,7 @@ Foam::radiation::greyMeanAbsorptionEmission::aCont(const label bandI) const
 Foam::tmp<Foam::volScalarField>
 Foam::radiation::greyMeanAbsorptionEmission::eCont(const label bandI) const
 {
-    tmp<volScalarField> e
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "eCont" + name(bandI),
-                mesh().time().timeName(),
-                mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh(),
-            dimensionedScalar("e", dimless/dimLength, 0.0)
-        )
-    );
-
-    return e;
+   return aCont(bandI);
 }
 
 
diff --git a/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C b/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C
index c2640f402a..dcaf7f01f6 100644
--- a/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C
+++ b/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C
@@ -231,24 +231,7 @@ Foam::radiation::wideBandAbsorptionEmission::aCont(const label bandI) const
 Foam::tmp<Foam::volScalarField>
 Foam::radiation::wideBandAbsorptionEmission::eCont(const label bandI) const
 {
-    tmp<volScalarField> e
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "e",
-                mesh().time().timeName(),
-                mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh(),
-            dimensionedScalar("e", dimless/dimLength, 0.0)
-        )
-    );
-
-    return e;
+    return aCont(bandI);
 }
 
 
diff --git a/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.C b/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.C
index 51120e89d7..e4032d854c 100644
--- a/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.C
+++ b/src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.C
@@ -268,7 +268,7 @@ kOmegaSST::kOmegaSST
       / max
         (
             a1_*omega_,
-            F2()*sqrt(2.0*magSqr(symm(fvc::grad(U_))))
+            F2()*sqrt(2.0)*mag(symm(fvc::grad(U_)))
         )
     );
     mut_.correctBoundaryConditions();
@@ -365,7 +365,7 @@ void kOmegaSST::correct()
         // Re-calculate viscosity
         mut_ =
             a1_*rho_*k_
-           /max(a1_*omega_, F2()*sqrt(2.0*magSqr(symm(fvc::grad(U_)))));
+           /max(a1_*omega_, F2()*sqrt(2.0)*mag(symm(fvc::grad(U_))));
         mut_.correctBoundaryConditions();
 
         // Re-calculate thermal diffusivity
diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/radiationProperties b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/radiationProperties
index 3282b1b075..ca98125036 100644
--- a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/radiationProperties
+++ b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/radiationProperties
@@ -137,6 +137,61 @@ greyMeanAbsorptionEmissionCoeffs
         );
     }
 
+    O2
+    {
+        Tcommon         300;
+        invTemp         true;
+        Tlow            200;
+        Thigh           2500;
+
+        loTcoeffs
+        (
+            0
+            0
+            0
+            0
+            0
+            0
+        );
+        hiTcoeffs
+        (
+            0.1
+            0
+            0
+            0
+            0
+            0
+        );
+    }
+
+
+    N2
+    {
+        Tcommon         300;
+        invTemp         true;
+        Tlow            200;
+        Thigh           2500;
+
+        loTcoeffs
+        (
+            0
+            0
+            0
+            0
+            0
+            0
+        );
+        hiTcoeffs
+        (
+            0.1
+            0
+            0
+            0
+            0
+            0
+        );
+    }
+
 }
 
 scatterModel    constantScatter;
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/radiationProperties b/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/radiationProperties
index 54fb40ca34..f9ff27414e 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/radiationProperties
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/radiationProperties
@@ -43,7 +43,7 @@ absorptionEmissionModel greyMeanAbsorptionEmission;
 constantAbsorptionEmissionCoeffs
 {
     a               a [ 0 -1 0 0 0 0 0 ] 0.01;
-    e               e [ 0 -1 0 0 0 0 0 ] 0;
+    e               e [ 0 -1 0 0 0 0 0 ] 0.01;
     E               E [ 1 -1 -3 0 0 0 0 ] 0;
 }
 
@@ -134,6 +134,61 @@ greyMeanAbsorptionEmissionCoeffs
             0
         );
     }
+
+    O2
+    {
+        Tcommon         300;
+        invTemp         true;
+        Tlow            200;
+        Thigh           2500;
+
+        loTcoeffs
+        (
+            0
+            0
+            0
+            0
+            0
+            0
+        );
+        hiTcoeffs
+        (
+            0.1
+            0
+            0
+            0
+            0
+            0
+        );
+    }
+
+
+    N2
+    {
+        Tcommon         300;
+        invTemp         true;
+        Tlow            200;
+        Thigh           2500;
+
+        loTcoeffs
+        (
+            0
+            0
+            0
+            0
+            0
+            0
+        );
+        hiTcoeffs
+        (
+            0.1
+            0
+            0
+            0
+            0
+            0
+        );
+    }
 }
 
 scatterModel    constantScatter;
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/radiationProperties b/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/radiationProperties
index 05fda138af..60bbcc81c4 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/radiationProperties
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/radiationProperties
@@ -135,6 +135,61 @@ greyMeanAbsorptionEmissionCoeffs
             0
         );
     }
+
+    O2
+    {
+        Tcommon         300;
+        invTemp         true;
+        Tlow            200;
+        Thigh           2500;
+
+        loTcoeffs
+        (
+            0
+            0
+            0
+            0
+            0
+            0
+        );
+        hiTcoeffs
+        (
+            0.1
+            0
+            0
+            0
+            0
+            0
+        );
+    }
+
+
+    N2
+    {
+        Tcommon         300;
+        invTemp         true;
+        Tlow            200;
+        Thigh           2500;
+
+        loTcoeffs
+        (
+            0
+            0
+            0
+            0
+            0
+            0
+        );
+        hiTcoeffs
+        (
+            0.1
+            0
+            0
+            0
+            0
+            0
+        );
+    }
 }
 
 scatterModel    constantScatter;
diff --git a/tutorials/mesh/cvMesh/blob/Allrun b/tutorials/mesh/cvMesh/blob/Allrun
index 1812d8e273..aafdcec484 100755
--- a/tutorials/mesh/cvMesh/blob/Allrun
+++ b/tutorials/mesh/cvMesh/blob/Allrun
@@ -6,7 +6,6 @@ cd ${0%/*} || exit 1    # run from this directory
 
 runApplication blockMesh
 runApplication cvMesh
-runApplication collapseEdges 1e-3 45
 runApplication snappyHexMesh
 runApplication checkMesh -constant -allGeometry -allTopology
 
diff --git a/tutorials/mesh/cvMesh/blob/system/cvMeshDict b/tutorials/mesh/cvMesh/blob/system/cvMeshDict
index b92aa466de..ca2fd3ccd6 100644
--- a/tutorials/mesh/cvMesh/blob/system/cvMeshDict
+++ b/tutorials/mesh/cvMesh/blob/system/cvMeshDict
@@ -59,7 +59,6 @@ initialPoints
 
     autoDensityCoeffs
     {
-        minCellSizeLimit 0.04;
         minLevels 0;
         maxSizeRatio 5.0;
         sampleResolution 5;
@@ -150,21 +149,20 @@ motionControl
         {
             priority            1;
             mode                bothSides;
+            
             surfaceCellSizeFunction uniformValue;
             uniformValueCoeffs
             {
                 surfaceCellSize 0.04;
             }
-            refinementFactor    1;
 
-            cellSizeFunction    uniformDistance;
+            cellSizeFunction        uniformDistance;
             uniformDistanceCoeffs
             {
                 distance        0.1;
             }
             linearDistanceCoeffs
             {
-                surfaceCellSize         0.02;
                 distanceCellSize        $defaultCellSize;
                 distance                0.1;
             }
@@ -217,14 +215,15 @@ polyMeshFiltering
 {
     writeTetDualMesh                        false;
     filterSizeCoeff                         0.2;
-    mergeClosenessCoeff                     1e-9;
+    mergeClosenessCoeff                     1e-3;
+    edgeMergeAngle                          45;
     continueFilteringOnBadInitialPolyMesh   true;
     filterErrorReductionCoeff               0.5;
     filterCountSkipThreshold                4;
     surfaceStepFaceAngle                    80;
     maxCollapseIterations                   25;
     maxConsecutiveEqualFaceSets             5;
-    edgeCollapseGuardFraction               0.3;
+    edgeCollapseGuardFraction               0.1;
     maxCollapseFaceToPointSideLengthCoeff   0.35;
 }
 
diff --git a/tutorials/mesh/cvMesh/blob/system/decomposeParDict b/tutorials/mesh/cvMesh/blob/system/decomposeParDict
index c347636c49..58cbdde811 100644
--- a/tutorials/mesh/cvMesh/blob/system/decomposeParDict
+++ b/tutorials/mesh/cvMesh/blob/system/decomposeParDict
@@ -17,8 +17,8 @@ FoamFile
 
 numberOfSubdomains 8;
 
-method          scotch;
-// method          ptscotch;
+//method          scotch;
+ method          ptscotch;
 // method          hierarchical;
 
 simpleCoeffs