diff --git a/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshTemplates.C b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshTemplates.C
index 63cef912b1..e5deecdf11 100644
--- a/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshTemplates.C
+++ b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshTemplates.C
@@ -65,8 +65,7 @@ void Foam::extendedFeatureEdgeMesh::sortPointsAndEdges
     labelListList featurePointFeatureEdges(nFeatPts);
     forAll(featurePointFeatureEdges, pI)
     {
-        featurePointFeatureEdges[pI] =
-            labelList(pointEdges[featurePoints[pI]].size(), -1);
+        featurePointFeatureEdges[pI] = pointEdges[featurePoints[pI]];
     }
 
     // Mapping between old and new indices, there is entry in the map for each
@@ -74,6 +73,10 @@ void Foam::extendedFeatureEdgeMesh::sortPointsAndEdges
     // >= 0 corresponds to the index
     labelList pointMap(sFeatLocalPts.size(), -1);
 
+    // Mapping between surface edge index and its feature edge index. -1 if it
+    // is not a feature edge
+    labelList edgeMap(sFeatEds.size(), -1);
+
     // Noting when the normal of a face has been used so not to duplicate
     labelList faceMap(surf.size(), -1);
 
@@ -98,6 +101,8 @@ void Foam::extendedFeatureEdgeMesh::sortPointsAndEdges
     {
         label sFEI = featureEdges[i];
 
+        edgeMap[sFEI] = i;
+
         const edge& fE(sFeatEds[sFEI]);
 
         // Check to see if the points have been already used
@@ -156,43 +161,31 @@ void Foam::extendedFeatureEdgeMesh::sortPointsAndEdges
         {
             regionEdges.append(i);
         }
-
-        forAll(featurePointFeatureEdges, pI)
-        {
-            const labelList& fpfEdges = pointEdges[featurePoints[pI]];
-
-            labelList& fpfe = featurePointFeatureEdges[pI];
-
-            forAll(fpfEdges, eI)
-            {
-                if (sFEI == fpfEdges[eI])
-                {
-                    fpfe[eI] = i;
-                }
-            }
-        }
     }
 
+    // Populate feature point feature edges
+    DynamicList<label> newFeatureEdges;
+
     forAll(featurePointFeatureEdges, pI)
     {
         const labelList& fpfe = featurePointFeatureEdges[pI];
 
-        DynamicList<label> newFeatureEdges(fpfe.size());
+        newFeatureEdges.setCapacity(fpfe.size());
 
         forAll(fpfe, eI)
         {
-            const label edgeIndex = fpfe[eI];
+            const label oldEdgeIndex = fpfe[eI];
+            const label newFeatureEdgeIndex = edgeMap[oldEdgeIndex];
 
-            if (edgeIndex != -1)
+            if (newFeatureEdgeIndex != -1)
             {
-                newFeatureEdges.append(edgeIndex);
+                newFeatureEdges.append(newFeatureEdgeIndex);
             }
         }
 
-        featurePointFeatureEdges[pI] = newFeatureEdges;
+        featurePointFeatureEdges[pI].transfer(newFeatureEdges);
     }
 
-
     // Reorder the edges by classification
 
     List<DynamicList<label> > allEds(nEdgeTypes);
@@ -277,7 +270,7 @@ void Foam::extendedFeatureEdgeMesh::sortPointsAndEdges
     edgeMesh::operator=(edgeMesh(pts, eds));
 
     // Initialise sorted edge related data
-    edgeDirections_ = edgeDirections/mag(edgeDirections);
+    edgeDirections_ = edgeDirections/(mag(edgeDirections) + VSMALL);
     edgeNormals_ = edgeNormals;
     regionEdges_ = regionEdges;
 
diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index 1ce1f78949..dd37136533 100644
--- a/src/finiteVolume/Make/files
+++ b/src/finiteVolume/Make/files
@@ -356,6 +356,8 @@ $(snGradSchemes)/faceCorrectedSnGrad/faceCorrectedSnGrads.C
 $(snGradSchemes)/limitedSnGrad/limitedSnGrads.C
 $(snGradSchemes)/uncorrectedSnGrad/uncorrectedSnGrads.C
 $(snGradSchemes)/orthogonalSnGrad/orthogonalSnGrads.C
+$(snGradSchemes)/quadraticFitSnGrad/quadraticFitSnGrads.C
+$(snGradSchemes)/linearFitSnGrad/linearFitSnGrads.C
 
 convectionSchemes = finiteVolume/convectionSchemes
 $(convectionSchemes)/convectionScheme/convectionSchemes.C
diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradData.C b/src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradData.C
new file mode 100644
index 0000000000..36a363ea1f
--- /dev/null
+++ b/src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradData.C
@@ -0,0 +1,270 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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 "CentredFitSnGradData.H"
+#include "surfaceFields.H"
+#include "SVD.H"
+#include "extendedCentredCellToFaceStencil.H"
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+template<class Polynomial>
+Foam::CentredFitSnGradData<Polynomial>::CentredFitSnGradData
+(
+    const fvMesh& mesh,
+    const extendedCentredCellToFaceStencil& stencil,
+    const scalar linearLimitFactor,
+    const scalar centralWeight
+)
+:
+    FitData
+    <
+        CentredFitSnGradData<Polynomial>,
+        extendedCentredCellToFaceStencil,
+        Polynomial
+    >
+    (
+        mesh, stencil, true, linearLimitFactor, centralWeight
+    ),
+    coeffs_(mesh.nFaces())
+{
+    if (debug)
+    {
+        Info<< "Contructing CentredFitSnGradData<Polynomial>" << endl;
+    }
+
+    calcFit();
+
+    if (debug)
+    {
+        Info<< "CentredFitSnGradData<Polynomial>::CentredFitSnGradData() :"
+            << "Finished constructing polynomialFit data"
+            << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Polynomial>
+void Foam::CentredFitSnGradData<Polynomial>::calcFit
+(
+    scalarList& coeffsi,
+    const List<point>& C,
+    const scalar wLin,
+    const scalar deltaCoeff,
+    const label facei
+)
+{
+    vector idir(1,0,0);
+    vector jdir(0,1,0);
+    vector kdir(0,0,1);
+    this->findFaceDirs(idir, jdir, kdir, facei);
+
+    // Setup the point weights
+    scalarList wts(C.size(), scalar(1));
+    wts[0] = this->centralWeight();
+    wts[1] = this->centralWeight();
+
+    // Reference point
+    point p0 = this->mesh().faceCentres()[facei];
+
+    // p0 -> p vector in the face-local coordinate system
+    vector d;
+
+    // Local coordinate scaling
+    scalar scale = 1;
+
+    // Matrix of the polynomial components
+    scalarRectangularMatrix B(C.size(), this->minSize(), scalar(0));
+
+    forAll(C, ip)
+    {
+        const point& p = C[ip];
+        const vector p0p = p - p0;
+
+        d.x() = p0p & idir;
+        d.y() = p0p & jdir;
+        d.z() = p0p & kdir;
+
+        if (ip == 0)
+        {
+            scale = cmptMax(cmptMag((d)));
+        }
+
+        // Scale the radius vector
+        d /= scale;
+
+        Polynomial::addCoeffs(B[ip], d, wts[ip], this->dim());
+    }
+
+    // Additional weighting for constant and linear terms
+    for (label i = 0; i < B.n(); i++)
+    {
+        B[i][0] *= wts[0];
+        B[i][1] *= wts[0];
+    }
+
+    // Set the fit
+    label stencilSize = C.size();
+    coeffsi.setSize(stencilSize);
+
+    bool goodFit = false;
+    for (int iIt = 0; iIt < 8 && !goodFit; iIt++)
+    {
+        SVD svd(B, SMALL);
+
+        for (label i=0; i<stencilSize; i++)
+        {
+            coeffsi[i] = wts[1]*wts[i]*svd.VSinvUt()[1][i]/scale;
+        }
+
+        goodFit =
+        (
+            mag(wts[0]*wts[0]*svd.VSinvUt()[0][0] - wLin)
+          < this->linearLimitFactor()*wLin)
+         && (mag(wts[0]*wts[1]*svd.VSinvUt()[0][1] - (1 - wLin)
+        ) < this->linearLimitFactor()*(1 - wLin))
+         && coeffsi[0] < 0 && coeffsi[1] > 0
+         && mag(coeffsi[0] + deltaCoeff) < 0.5*deltaCoeff
+         && mag(coeffsi[1] - deltaCoeff) < 0.5*deltaCoeff;
+
+        if (!goodFit)
+        {
+            // (not good fit so increase weight in the centre and weight
+            //  for constant and linear terms)
+
+            WarningIn
+            (
+                "CentredFitSnGradData<Polynomial>::calcFit"
+                "(const List<point>& C, const label facei"
+            )   << "Cannot fit face " << facei << " iteration " << iIt
+                << " with sum of weights " << sum(coeffsi) << nl
+                << "    Weights " << coeffsi << nl
+                << "    Linear weights " << wLin << " " << 1 - wLin << nl
+                << "    deltaCoeff " << deltaCoeff << nl
+                << "    sing vals " << svd.S() << nl
+                << "Components of goodFit:\n"
+                << "    wts[0]*wts[0]*svd.VSinvUt()[0][0] = "
+                << wts[0]*wts[0]*svd.VSinvUt()[0][0] << nl
+                << "    wts[0]*wts[1]*svd.VSinvUt()[0][1] = "
+                << wts[0]*wts[1]*svd.VSinvUt()[0][1]
+                << " dim = " << this->dim() << endl;
+
+            wts[0] *= 10;
+            wts[1] *= 10;
+
+            for (label j = 0; j < B.m(); j++)
+            {
+                B[0][j] *= 10;
+                B[1][j] *= 10;
+            }
+
+            for (label i = 0; i < B.n(); i++)
+            {
+                B[i][0] *= 10;
+                B[i][1] *= 10;
+            }
+        }
+    }
+
+    if (goodFit)
+    {
+        // Remove the uncorrected coefficients
+        coeffsi[0] += deltaCoeff;
+        coeffsi[1] -= deltaCoeff;
+    }
+    else
+    {
+        WarningIn
+        (
+            "CentredFitSnGradData<Polynomial>::calcFit(..)"
+        )   << "Could not fit face " << facei
+            << "    Coefficients = " << coeffsi
+            << ", reverting to uncorrected." << endl;
+
+        coeffsi = 0;
+    }
+}
+
+
+template<class Polynomial>
+void Foam::CentredFitSnGradData<Polynomial>::calcFit()
+{
+    const fvMesh& mesh = this->mesh();
+
+    // Get the cell/face centres in stencil order.
+    // Centred face stencils no good for triangles or tets.
+    // Need bigger stencils
+    List<List<point> > stencilPoints(mesh.nFaces());
+    this->stencil().collectData(mesh.C(), stencilPoints);
+
+    // find the fit coefficients for every face in the mesh
+
+    const surfaceScalarField& w = mesh.surfaceInterpolation::weights();
+    const surfaceScalarField& dC = mesh.nonOrthDeltaCoeffs();
+
+    for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
+    {
+        calcFit
+        (
+            coeffs_[facei],
+            stencilPoints[facei],
+            w[facei],
+            dC[facei],
+            facei
+        );
+    }
+
+    const surfaceScalarField::GeometricBoundaryField& bw = w.boundaryField();
+    const surfaceScalarField::GeometricBoundaryField& bdC = dC.boundaryField();
+
+    forAll(bw, patchi)
+    {
+        const fvsPatchScalarField& pw = bw[patchi];
+        const fvsPatchScalarField& pdC = bdC[patchi];
+
+        if (pw.coupled())
+        {
+            label facei = pw.patch().start();
+
+            forAll(pw, i)
+            {
+                calcFit
+                (
+                    coeffs_[facei],
+                    stencilPoints[facei],
+                    pw[i],
+                    pdC[i],
+                    facei
+                );
+                facei++;
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradData.H b/src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradData.H
new file mode 100644
index 0000000000..be834684a7
--- /dev/null
+++ b/src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradData.H
@@ -0,0 +1,127 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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::CentredFitSnGradData
+
+Description
+    Data for centred fit snGrad schemes
+
+SourceFiles
+    CentredFitSnGradData.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef CentredFitSnGradData_H
+#define CentredFitSnGradData_H
+
+#include "FitData.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class extendedCentredCellToFaceStencil;
+
+/*---------------------------------------------------------------------------*\
+                    Class CentredFitSnGradData Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Polynomial>
+class CentredFitSnGradData
+:
+    public FitData
+    <
+        CentredFitSnGradData<Polynomial>,
+        extendedCentredCellToFaceStencil,
+        Polynomial
+    >
+{
+    // Private data
+
+        //- For each cell in the mesh store the values which multiply the
+        //  values of the stencil to obtain the gradient for each direction
+        List<scalarList> coeffs_;
+
+
+public:
+
+    TypeName("CentredFitSnGradData");
+
+
+    // Constructors
+
+        //- Construct from components
+        CentredFitSnGradData
+        (
+            const fvMesh& mesh,
+            const extendedCentredCellToFaceStencil& stencil,
+            const scalar linearLimitFactor,
+            const scalar centralWeight
+        );
+
+
+    //- Destructor
+    virtual ~CentredFitSnGradData()
+    {}
+
+
+    // Member functions
+
+        //- Return reference to fit coefficients
+        const List<scalarList>& coeffs() const
+        {
+            return coeffs_;
+        }
+
+        //- Calculate the fit for the specified face and set the coefficients
+        void calcFit
+        (
+            scalarList& coeffsi, // coefficients to be set
+            const List<point>&,  // Stencil points
+            const scalar wLin,   // Weight for linear approximation (weights
+                                 // nearest neighbours)
+            const scalar deltaCoeff, // uncorrected delta coefficient
+            const label faci     // Current face index
+        );
+
+        void calcFit();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "CentredFitSnGradData.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradScheme.H b/src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradScheme.H
new file mode 100644
index 0000000000..d92478c317
--- /dev/null
+++ b/src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradScheme.H
@@ -0,0 +1,183 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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::CentredFitSnGradScheme
+
+Description
+    Centred fit snGrad scheme which applies an explicit correction to snGrad
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef CentredFitSnGradScheme_H
+#define CentredFitSnGradScheme_H
+
+#include "CentredFitSnGradData.H"
+#include "snGradScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class fvMesh;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fv
+{
+/*---------------------------------------------------------------------------*\
+                           Class CentredFitSnGradSnGradScheme Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type, class Polynomial, class Stencil>
+class CentredFitSnGradScheme
+:
+    public snGradScheme<Type>
+{
+    // Private Data
+
+        //- Factor the fit is allowed to deviate from linear.
+        //  This limits the amount of high-order correction and increases
+        //  stability on bad meshes
+        const scalar linearLimitFactor_;
+
+        //- Weights for central stencil
+        const scalar centralWeight_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        CentredFitSnGradScheme(const CentredFitSnGradScheme&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const CentredFitSnGradScheme&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("CentredFitSnGradScheme");
+
+
+    // Constructors
+
+        //- Construct from mesh and Istream
+        CentredFitSnGradScheme(const fvMesh& mesh, Istream& is)
+        :
+            snGradScheme<Type>(mesh),
+            linearLimitFactor_(readScalar(is)),
+            centralWeight_(1000)
+        {}
+
+
+    // Member Functions
+
+        //- Return the interpolation weighting factors for the given field
+        virtual tmp<surfaceScalarField> deltaCoeffs
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const
+        {
+            return this->mesh().nonOrthDeltaCoeffs();
+        }
+
+        //- Return true if this scheme uses an explicit correction
+        virtual bool corrected() const
+        {
+            return true;
+        }
+
+        //- Return the explicit correction to the face-interpolate
+        virtual tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
+        correction
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& vf
+        ) const
+        {
+            const fvMesh& mesh = this->mesh();
+
+            const extendedCentredCellToFaceStencil& stencil = Stencil::New
+            (
+                mesh
+            );
+
+            const CentredFitSnGradData<Polynomial>& cfd =
+                CentredFitSnGradData<Polynomial>::New
+                (
+                    mesh,
+                    stencil,
+                    linearLimitFactor_,
+                    centralWeight_
+                );
+
+            tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > sft
+            (
+                stencil.weightedSum(vf, cfd.coeffs())
+            );
+
+            sft().dimensions() /= dimLength;
+
+            return sft;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Add the patch constructor functions to the hash tables
+
+#define makeCentredFitSnGradTypeScheme(SS, POLYNOMIAL, STENCIL, TYPE) \
+                                                                      \
+typedef CentredFitSnGradScheme<TYPE, POLYNOMIAL, STENCIL>             \
+    CentredFitSnGradScheme##TYPE##POLYNOMIAL##STENCIL##_;             \
+defineTemplateTypeNameAndDebugWithName                                \
+    (CentredFitSnGradScheme##TYPE##POLYNOMIAL##STENCIL##_, #SS, 0);   \
+                                                                      \
+snGradScheme<TYPE>::addMeshConstructorToTable                         \
+<CentredFitSnGradScheme<TYPE, POLYNOMIAL, STENCIL> >                  \
+    add##SS##STENCIL##TYPE##MeshConstructorToTable_;
+
+#define makeCentredFitSnGradScheme(SS, POLYNOMIAL, STENCIL)           \
+                                                                      \
+makeCentredFitSnGradTypeScheme(SS,POLYNOMIAL,STENCIL,scalar)          \
+makeCentredFitSnGradTypeScheme(SS,POLYNOMIAL,STENCIL,vector)          \
+makeCentredFitSnGradTypeScheme(SS,POLYNOMIAL,STENCIL,sphericalTensor) \
+makeCentredFitSnGradTypeScheme(SS,POLYNOMIAL,STENCIL,symmTensor)      \
+makeCentredFitSnGradTypeScheme(SS,POLYNOMIAL,STENCIL,tensor)
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/linearFitSnGrad/linearFitSnGrads.C b/src/finiteVolume/finiteVolume/snGradSchemes/linearFitSnGrad/linearFitSnGrads.C
new file mode 100644
index 0000000000..4335cc8e69
--- /dev/null
+++ b/src/finiteVolume/finiteVolume/snGradSchemes/linearFitSnGrad/linearFitSnGrads.C
@@ -0,0 +1,51 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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 "CentredFitSnGradScheme.H"
+#include "linearFitPolynomial.H"
+#include "centredFECCellToFaceStencilObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTemplateTypeNameAndDebug
+    (
+        CentredFitSnGradData<linearFitPolynomial>,
+        0
+    );
+
+    namespace fv
+    {
+        makeCentredFitSnGradScheme
+        (
+            linearFit,
+            linearFitPolynomial,
+            centredFECCellToFaceStencilObject
+        );
+    }
+}
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGrads.C b/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGrads.C
new file mode 100644
index 0000000000..63e4ea98fa
--- /dev/null
+++ b/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGrads.C
@@ -0,0 +1,51 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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 "CentredFitSnGradScheme.H"
+#include "quadraticFitPolynomial.H"
+#include "centredCFCCellToFaceStencilObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTemplateTypeNameAndDebug
+    (
+        CentredFitSnGradData<quadraticFitPolynomial>,
+        0
+    );
+
+    namespace fv
+    {
+        makeCentredFitSnGradScheme
+        (
+            quadraticFit,
+            quadraticFitPolynomial,
+            centredCFCCellToFaceStencilObject
+        );
+    }
+}
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C
index 5b227e12eb..c1d322ff4e 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C
@@ -45,11 +45,11 @@ Foam::FitData<Form, ExtendedStencil, Polynomial>::FitData
     linearCorrection_(linearCorrection),
     linearLimitFactor_(linearLimitFactor),
     centralWeight_(centralWeight),
-#   ifdef SPHERICAL_GEOMETRY
+    #ifdef SPHERICAL_GEOMETRY
     dim_(2),
-#   else
+    #else
     dim_(mesh.nGeometricD()),
-#   endif
+    #endif
     minSize_(Polynomial::nTerms(dim_))
 {
     // Check input
@@ -79,7 +79,7 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::findFaceDirs
     idir = mesh.faceAreas()[facei];
     idir /= mag(idir);
 
-#   ifndef SPHERICAL_GEOMETRY
+    #ifndef SPHERICAL_GEOMETRY
     if (mesh.nGeometricD() <= 2) // find the normal direction
     {
         if (mesh.geometricD()[0] == -1)
@@ -100,10 +100,10 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::findFaceDirs
         const face& f = mesh.faces()[facei];
         kdir = mesh.points()[f[0]] - mesh.faceCentres()[facei];
     }
-#   else
+    #else
     // Spherical geometry so kdir is the radial direction
     kdir = mesh.faceCentres()[facei];
-#   endif
+    #endif
 
     if (mesh.nGeometricD() == 3)
     {
@@ -170,11 +170,11 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::calcFit
 
         d.x() = (p - p0)&idir;
         d.y() = (p - p0)&jdir;
-#       ifndef SPHERICAL_GEOMETRY
+        #ifndef SPHERICAL_GEOMETRY
         d.z() = (p - p0)&kdir;
-#       else
+        #else
         d.z() = mag(p) - mag(p0);
-#       endif
+        #endif
 
         if (ip == 0)
         {
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.H
index cee373b6c4..5faeed6eb1 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.H
@@ -29,7 +29,7 @@ Description
     The linearCorrection_ determines whether the fit is for a corrected
     linear scheme (first two coefficients are corrections for owner and
     neighbour) or a pure upwind scheme (first coefficient is correction for
-    owner ; weight on face taken as 1).
+    owner; weight on face taken as 1).
 
 SourceFiles
     FitData.C
@@ -80,7 +80,7 @@ class FitData
         const label minSize_;
 
 
-    // Private Member Functions
+protected:
 
         //- Find the normal direction (i) and j and k directions for face faci
         void findFaceDirs
@@ -93,9 +93,6 @@ class FitData
 
 public:
 
-    //TypeName("FitData");
-
-
     // Constructors
 
         //- Construct from components
@@ -122,6 +119,30 @@ public:
             return stencil_;
         }
 
+        //- Factor the fit is allowed to deviate from the base scheme
+        scalar linearLimitFactor() const
+        {
+            return linearLimitFactor_;
+        }
+
+        //- Return weight for central stencil
+        scalar centralWeight() const
+        {
+            return centralWeight_;
+        }
+
+        //- Dimensionality of the geometry
+        label dim() const
+        {
+            return dim_;
+        }
+
+        //- Minimum stencil size
+        label minSize() const
+        {
+            return minSize_;
+        }
+
         bool linearCorrection() const
         {
             return linearCorrection_;
@@ -140,7 +161,6 @@ public:
         //- Calculate the fit for all the faces
         virtual void calcFit() = 0;
 
-
         //- Recalculate weights (but not stencil) when the mesh moves
         bool movePoints();
 };
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
index cecee33cef..37e61cbbb7 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
@@ -767,7 +767,7 @@ void Foam::autoSnapDriver::binFeatureFace
 ) const
 {
     // What to do with very far attraction? For now just ignore the face
-    if (magSqr(faceDisp) < sqr(snapDist))
+    if (magSqr(faceDisp) < sqr(snapDist) && mag(faceSurfaceNormal) > VSMALL)
     {
         const point pt = fc + faceDisp;
 
diff --git a/src/thermophysicalModels/radiationModels/radiationModel/P1/P1.C b/src/thermophysicalModels/radiationModels/radiationModel/P1/P1.C
index a07142897b..cf7a2329d4 100644
--- a/src/thermophysicalModels/radiationModels/radiationModel/P1/P1.C
+++ b/src/thermophysicalModels/radiationModels/radiationModel/P1/P1.C
@@ -241,10 +241,14 @@ void Foam::radiation::P1::calculate()
     );
 
     // Calculate radiative heat flux on boundaries.
-    forAll(mesh_.boundaryMesh(), patchI)
+    forAll(mesh_.boundaryMesh(), patchi)
     {
-        Qr_.boundaryField()[patchI] =
-            -gamma.boundaryField()[patchI]*G_.boundaryField()[patchI].snGrad();
+        if (!G_.boundaryField()[patchi].coupled())
+        {
+            Qr_.boundaryField()[patchi] =
+                -gamma.boundaryField()[patchi]
+                *G_.boundaryField()[patchi].snGrad();
+        }
     }
 }
 
diff --git a/src/transportModels/incompressible/incompressibleTwoPhaseMixture/incompressibleTwoPhaseMixture.C b/src/transportModels/incompressible/incompressibleTwoPhaseMixture/incompressibleTwoPhaseMixture.C
index 98046c9b76..397f9566b1 100644
--- a/src/transportModels/incompressible/incompressibleTwoPhaseMixture/incompressibleTwoPhaseMixture.C
+++ b/src/transportModels/incompressible/incompressibleTwoPhaseMixture/incompressibleTwoPhaseMixture.C
@@ -175,8 +175,14 @@ bool Foam::incompressibleTwoPhaseMixture::read()
     {
         if
         (
-            nuModel1_().read(subDict(phase1Name_))
-         && nuModel2_().read(subDict(phase2Name_))
+            nuModel1_().read
+            (
+                subDict(phase1Name_ == "1" ? "phase1": phase1Name_)
+            )
+         && nuModel2_().read
+            (
+                subDict(phase2Name_ == "2" ? "phase2": phase2Name_)
+            )
         )
         {
             nuModel1_->viscosityProperties().lookup("rho") >> rho1_;