diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C
index 0191a90a8c..d9898beda1 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C
@@ -55,6 +55,64 @@ scalar epsilonLowReWallFunctionFvPatchScalarField::yPlusLam
 }
 
 
+void epsilonLowReWallFunctionFvPatchScalarField::calculate
+(
+    const turbulenceModel& turbulence,
+    const List<scalar>& cornerWeights,
+    const fvPatch& patch,
+    scalarField& G,
+    scalarField& epsilon
+)
+{
+    const label patchI = patch.index();
+
+    const scalarField& y = turbulence.y()[patchI];
+
+    const scalar Cmu25 = pow025(Cmu_);
+    const scalar Cmu75 = pow(Cmu_, 0.75);
+
+    const tmp<volScalarField> tk = turbulence.k();
+    const volScalarField& k = tk();
+
+    const tmp<volScalarField> tnu = turbulence.nu();
+    const scalarField& nuw = tnu().boundaryField()[patchI];
+
+    const tmp<volScalarField> tnut = turbulence.nut();
+    const volScalarField& nut = tnut();
+    const scalarField& nutw = nut.boundaryField()[patchI];
+
+    const fvPatchVectorField& Uw = turbulence.U().boundaryField()[patchI];
+
+    const scalarField magGradUw(mag(Uw.snGrad()));
+
+    // Set epsilon and G
+    forAll(nutw, faceI)
+    {
+        label cellI = patch.faceCells()[faceI];
+
+        scalar yPlus = Cmu25*sqrt(k[cellI])*y[faceI]/nuw[faceI];
+
+        scalar w = cornerWeights[faceI];
+
+        if (yPlus > yPlusLam_)
+        {
+            epsilon[cellI] = w*Cmu75*pow(k[cellI], 1.5)/(kappa_*y[faceI]);
+        }
+        else
+        {
+            epsilon[cellI] = w*2.0*k[cellI]*nuw[faceI]/sqr(y[faceI]);
+        }
+
+        G[cellI] =
+            w
+           *(nutw[faceI] + nuw[faceI])
+           *magGradUw[faceI]
+           *Cmu25*sqrt(k[cellI])
+           /(kappa_*y[faceI]);
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 epsilonLowReWallFunctionFvPatchScalarField::
@@ -119,81 +177,6 @@ epsilonLowReWallFunctionFvPatchScalarField
 {}
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void epsilonLowReWallFunctionFvPatchScalarField::updateCoeffs()
-{
-    if (updated())
-    {
-        return;
-    }
-
-    const label patchI = patch().index();
-
-    const turbulenceModel& turbulence =
-        db().lookupObject<turbulenceModel>("turbulenceModel");
-    const scalarField& y = turbulence.y()[patchI];
-
-    volScalarField& G =
-        const_cast<volScalarField&>
-        (
-            db().lookupObject<volScalarField>
-            (
-                turbulence.GName()
-            )
-        );
-
-    DimensionedField<scalar, volMesh>& epsilon =
-        const_cast<DimensionedField<scalar, volMesh>&>
-        (
-            dimensionedInternalField()
-        );
-
-    const tmp<volScalarField> tk = turbulence.k();
-    const volScalarField& k = tk();
-
-    const tmp<volScalarField> tnu = turbulence.nu();
-    const scalarField& nuw = tnu().boundaryField()[patchI];
-
-    const tmp<volScalarField> tnut = turbulence.nut();
-    const volScalarField& nut = tnut();
-    const scalarField& nutw = nut.boundaryField()[patchI];
-
-    const fvPatchVectorField& Uw = turbulence.U().boundaryField()[patchI];
-    const scalarField magGradUw(mag(Uw.snGrad()));
-
-    const scalar Cmu25 = pow025(Cmu_);
-    const scalar Cmu75 = pow(Cmu_, 0.75);
-
-    // Set epsilon and G
-    forAll(nutw, faceI)
-    {
-        label faceCellI = patch().faceCells()[faceI];
-
-        scalar yPlus = Cmu25*sqrt(k[faceCellI])*y[faceI]/nuw[faceI];
-
-        if (yPlus > yPlusLam_)
-        {
-            epsilon[faceCellI] = Cmu75*pow(k[faceCellI], 1.5)/(kappa_*y[faceI]);
-        }
-        else
-        {
-            epsilon[faceCellI] = 2.0*k[faceCellI]*nuw[faceI]/sqr(y[faceI]);
-        }
-
-        G[faceCellI] =
-            (nutw[faceI] + nuw[faceI])
-           *magGradUw[faceI]
-           *Cmu25*sqrt(k[faceCellI])
-           /(kappa_*y[faceI]);
-    }
-
-    fixedInternalValueFvPatchField<scalar>::updateCoeffs();
-
-    // TODO: perform averaging for cells sharing more than one boundary face
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 makePatchTypeField
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.H b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.H
index b8cf09f0b7..778fffcfd7 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.H
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -97,6 +97,16 @@ protected:
         //- Calculate the Y+ at the edge of the laminar sublayer
         scalar yPlusLam(const scalar kappa, const scalar E);
 
+        //- Calculate the epsilon and G
+        virtual void calculate
+        (
+            const turbulenceModel& turbulence,
+            const List<scalar>& cornerWeights,
+            const fvPatch& patch,
+            scalarField& G,
+            scalarField& epsilon
+        );
+
 
 public:
 
@@ -166,13 +176,9 @@ public:
             );
         }
 
-
-    // Member functions
-
-        // Evaluation functions
-
-            //- Update the coefficients associated with the patch field
-            virtual void updateCoeffs();
+    //- Destructor
+    virtual ~epsilonLowReWallFunctionFvPatchScalarField()
+    {}
 };
 
 
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C
index 9435962e6a..7fb6fcfb96 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C
@@ -26,11 +26,10 @@ License
 #include "epsilonWallFunctionFvPatchScalarField.H"
 #include "incompressible/turbulenceModel/turbulenceModel.H"
 #include "fvPatchFieldMapper.H"
+#include "fvMatrix.H"
 #include "volFields.H"
-#include "addToRunTimeSelectionTable.H"
 #include "wallFvPatch.H"
-#include "nutkWallFunctionFvPatchScalarField.H"
-
+#include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -39,7 +38,7 @@ namespace Foam
 namespace incompressible
 {
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
 
 void epsilonWallFunctionFvPatchScalarField::checkType()
 {
@@ -63,118 +62,160 @@ void epsilonWallFunctionFvPatchScalarField::writeLocalEntries(Ostream& os) const
 }
 
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedInternalValueFvPatchField<scalar>(p, iF),
-    Cmu_(0.09),
-    kappa_(0.41),
-    E_(9.8)
+void epsilonWallFunctionFvPatchScalarField::setMaster()
 {
-    checkType();
-}
-
-
-epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
-(
-    const epsilonWallFunctionFvPatchScalarField& ptf,
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const fvPatchFieldMapper& mapper
-)
-:
-    fixedInternalValueFvPatchField<scalar>(ptf, p, iF, mapper),
-    Cmu_(ptf.Cmu_),
-    kappa_(ptf.kappa_),
-    E_(ptf.E_)
-{
-    checkType();
-}
-
-
-epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const dictionary& dict
-)
-:
-    fixedInternalValueFvPatchField<scalar>(p, iF, dict),
-    Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
-    kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
-    E_(dict.lookupOrDefault<scalar>("E", 9.8))
-{
-    checkType();
-}
-
-
-epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
-(
-    const epsilonWallFunctionFvPatchScalarField& ewfpsf
-)
-:
-    fixedInternalValueFvPatchField<scalar>(ewfpsf),
-    Cmu_(ewfpsf.Cmu_),
-    kappa_(ewfpsf.kappa_),
-    E_(ewfpsf.E_)
-{
-    checkType();
-}
-
-
-epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
-(
-    const epsilonWallFunctionFvPatchScalarField& ewfpsf,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedInternalValueFvPatchField<scalar>(ewfpsf, iF),
-    Cmu_(ewfpsf.Cmu_),
-    kappa_(ewfpsf.kappa_),
-    E_(ewfpsf.E_)
-{
-    checkType();
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void epsilonWallFunctionFvPatchScalarField::updateCoeffs()
-{
-    if (updated())
+    if (master_ != -1)
     {
         return;
     }
 
-    const label patchI = patch().index();
+    const volScalarField& epsilon =
+        static_cast<const volScalarField&>(this->dimensionedInternalField());
+
+    const volScalarField::GeometricBoundaryField& bf = epsilon.boundaryField();
+
+    label master = -1;
+    forAll(bf, patchI)
+    {
+        if (isA<epsilonWallFunctionFvPatchScalarField>(bf[patchI]))
+        {
+            epsilonWallFunctionFvPatchScalarField& epf = epsilonPatch(patchI);
+
+            if (master == -1)
+            {
+                master = patchI;
+            }
+
+            epf.master() = master;
+        }
+    }
+}
+
+
+void epsilonWallFunctionFvPatchScalarField::createAveragingWeights()
+{
+    if (initialised_)
+    {
+        return;
+    }
+
+    const volScalarField& epsilon =
+        static_cast<const volScalarField&>(this->dimensionedInternalField());
+
+    const volScalarField::GeometricBoundaryField& bf = epsilon.boundaryField();
+
+    const fvMesh& mesh = epsilon.mesh();
+
+    volScalarField weights
+    (
+        IOobject
+        (
+            "weights",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false // do not register
+        ),
+        mesh,
+        dimensionedScalar("zero", dimless, 0.0)
+    );
+
+    DynamicList<label> epsilonPatches(bf.size());
+    forAll(bf, patchI)
+    {
+        if (isA<epsilonWallFunctionFvPatchScalarField>(bf[patchI]))
+        {
+            epsilonPatches.append(patchI);
+
+            const labelUList& faceCells = bf[patchI].patch().faceCells();
+            forAll(faceCells, i)
+            {
+                label cellI = faceCells[i];
+                weights[cellI]++;
+            }
+        }
+    }
+
+    cornerWeights_.setSize(bf.size());
+    forAll(epsilonPatches, i)
+    {
+        label patchI = epsilonPatches[i];
+        const fvPatchField& wf = weights.boundaryField()[patchI];
+        cornerWeights_[patchI] = 1.0/wf.patchInternalField();
+    }
+
+    G_.setSize(dimensionedInternalField().size(), 0.0);
+    epsilon_.setSize(dimensionedInternalField().size(), 0.0);
+
+    initialised_ = true;
+}
+
+
+epsilonWallFunctionFvPatchScalarField&
+epsilonWallFunctionFvPatchScalarField::epsilonPatch(const label patchI)
+{
+    const volScalarField& epsilon =
+        static_cast<const volScalarField&>(this->dimensionedInternalField());
+
+    const volScalarField::GeometricBoundaryField& bf = epsilon.boundaryField();
+
+    const epsilonWallFunctionFvPatchScalarField& epf =
+        refCast<const epsilonWallFunctionFvPatchScalarField>(bf[patchI]);
+
+    return const_cast<epsilonWallFunctionFvPatchScalarField&>(epf);
+}
+
+
+void epsilonWallFunctionFvPatchScalarField::calculateTurbulenceFields
+(
+    const turbulenceModel& turbulence,
+    scalarField& G0,
+    scalarField& epsilon0
+)
+{
+    // accumulate all of the G and epsilon contributions
+    forAll(cornerWeights_, patchI)
+    {
+        if (!cornerWeights_[patchI].empty())
+        {
+            epsilonWallFunctionFvPatchScalarField& epf = epsilonPatch(patchI);
+
+            const List<scalar>& w = cornerWeights_[patchI];
+
+            epf.calculate(turbulence, w, epf.patch(), G0, epsilon0);
+        }
+    }
+
+    // apply zero-gradient condition for epsilon
+    forAll(cornerWeights_, patchI)
+    {
+        if (!cornerWeights_[patchI].empty())
+        {
+            epsilonWallFunctionFvPatchScalarField& epf = epsilonPatch(patchI);
+
+            epf == scalarField(epsilon0, epf.patch().faceCells());
+        }
+    }
+}
+
+
+void epsilonWallFunctionFvPatchScalarField::calculate
+(
+    const turbulenceModel& turbulence,
+    const List<scalar>& cornerWeights,
+    const fvPatch& patch,
+    scalarField& G,
+    scalarField& epsilon
+)
+{
+    const label patchI = patch.index();
 
-    const turbulenceModel& turbulence =
-        db().lookupObject<turbulenceModel>("turbulenceModel");
     const scalarField& y = turbulence.y()[patchI];
 
     const scalar Cmu25 = pow025(Cmu_);
     const scalar Cmu75 = pow(Cmu_, 0.75);
 
-    volScalarField& G =
-        const_cast<volScalarField&>
-        (
-            db().lookupObject<volScalarField>
-            (
-                turbulence.GName()
-            )
-        );
-
-    DimensionedField<scalar, volMesh>& epsilon =
-        const_cast<DimensionedField<scalar, volMesh>&>
-        (
-            dimensionedInternalField()
-        );
-
     const tmp<volScalarField> tk = turbulence.k();
     const volScalarField& k = tk();
 
@@ -192,35 +233,329 @@ void epsilonWallFunctionFvPatchScalarField::updateCoeffs()
     // Set epsilon and G
     forAll(nutw, faceI)
     {
-        label faceCellI = patch().faceCells()[faceI];
+        label cellI = patch.faceCells()[faceI];
 
-        epsilon[faceCellI] = Cmu75*pow(k[faceCellI], 1.5)/(kappa_*y[faceI]);
+        scalar w = cornerWeights[faceI];
 
-        G[faceCellI] =
-            (nutw[faceI] + nuw[faceI])
+        epsilon[cellI] += w*Cmu75*pow(k[cellI], 1.5)/(kappa_*y[faceI]);
+
+        G[cellI] +=
+            w
+           *(nutw[faceI] + nuw[faceI])
            *magGradUw[faceI]
-           *Cmu25*sqrt(k[faceCellI])
+           *Cmu25*sqrt(k[cellI])
            /(kappa_*y[faceI]);
     }
-
-    fixedInternalValueFvPatchField<scalar>::updateCoeffs();
-
-    // TODO: perform averaging for cells sharing more than one boundary face
 }
 
 
-void epsilonWallFunctionFvPatchScalarField::evaluate
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
 (
-    const Pstream::commsTypes commsType
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchField<scalar>(p, iF),
+    Cmu_(0.09),
+    kappa_(0.41),
+    E_(9.8),
+    G_(),
+    epsilon_(),
+    initialised_(false),
+    master_(-1),
+    cornerWeights_()
+{
+    checkType();
+}
+
+
+epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
+(
+    const epsilonWallFunctionFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchField<scalar>(ptf, p, iF, mapper),
+    Cmu_(ptf.Cmu_),
+    kappa_(ptf.kappa_),
+    E_(ptf.E_),
+    G_(),
+    epsilon_(),
+    initialised_(false),
+    master_(-1),
+    cornerWeights_()
+{
+    checkType();
+}
+
+
+epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchField<scalar>(p, iF, dict),
+    Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
+    kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
+    E_(dict.lookupOrDefault<scalar>("E", 9.8)),
+    G_(),
+    epsilon_(),
+    initialised_(false),
+    master_(-1),
+    cornerWeights_()
+{
+    checkType();
+}
+
+
+epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
+(
+    const epsilonWallFunctionFvPatchScalarField& ewfpsf
+)
+:
+    fixedValueFvPatchField<scalar>(ewfpsf),
+    Cmu_(ewfpsf.Cmu_),
+    kappa_(ewfpsf.kappa_),
+    E_(ewfpsf.E_),
+    G_(),
+    epsilon_(),
+    initialised_(false),
+    master_(-1),
+    cornerWeights_()
+{
+    checkType();
+}
+
+
+epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
+(
+    const epsilonWallFunctionFvPatchScalarField& ewfpsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchField<scalar>(ewfpsf, iF),
+    Cmu_(ewfpsf.Cmu_),
+    kappa_(ewfpsf.kappa_),
+    E_(ewfpsf.E_),
+    G_(),
+    epsilon_(),
+    initialised_(false),
+    master_(-1),
+    cornerWeights_()
+{
+    checkType();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+scalarField& epsilonWallFunctionFvPatchScalarField::G(bool init)
+{
+    if (patch().index() == master_)
+    {
+        if (init)
+        {
+            G_ = 0.0;
+        }
+
+        return G_;
+    }
+
+    return epsilonPatch(master_).G();
+}
+
+
+scalarField& epsilonWallFunctionFvPatchScalarField::epsilon(bool init)
+{
+    if (patch().index() == master_)
+    {
+        if (init)
+        {
+            epsilon_ = 0.0;
+        }
+
+        return epsilon_;
+    }
+
+    return epsilonPatch(master_).epsilon(init);
+}
+
+
+void epsilonWallFunctionFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const turbulenceModel& turbulence =
+        db().lookupObject<turbulenceModel>(turbulenceModel::typeName);
+
+    setMaster();
+
+    if (patch().index() == master_)
+    {
+        createAveragingWeights();
+        calculateTurbulenceFields(turbulence, G(true), epsilon(true));
+    }
+
+    const scalarField& G0 = this->G();
+    const scalarField& epsilon0 = this->epsilon();
+
+    typedef DimensionedField<scalar, volMesh> FieldType;
+
+    FieldType& G =
+        const_cast<FieldType&>
+        (
+            db().lookupObject<FieldType>(turbulence.GName())
+        );
+
+    FieldType& epsilon = const_cast<FieldType&>(dimensionedInternalField());
+
+    forAll(*this, faceI)
+    {
+        label cellI = patch().faceCells()[faceI];
+
+        G[cellI] = G0[cellI];
+        epsilon[cellI] = epsilon0[cellI];
+    }
+
+    fvPatchField<scalar>::updateCoeffs();
+}
+
+
+void epsilonWallFunctionFvPatchScalarField::updateCoeffs
+(
+    const scalarField& weights
 )
 {
-    fixedInternalValueFvPatchField<scalar>::evaluate(commsType);
+    if (updated())
+    {
+        return;
+    }
+
+    const turbulenceModel& turbulence =
+        db().lookupObject<turbulenceModel>(turbulenceModel::typeName);
+
+    setMaster();
+
+    if (patch().index() == master_)
+    {
+        createAveragingWeights();
+        calculateTurbulenceFields(turbulence, G(true), epsilon(true));
+    }
+
+    const scalarField& G0 = this->G();
+    const scalarField& epsilon0 = this->epsilon();
+
+    typedef DimensionedField<scalar, volMesh> FieldType;
+
+    FieldType& G =
+        const_cast<FieldType&>
+        (
+            db().lookupObject<FieldType>(turbulence.GName())
+        );
+
+    FieldType& epsilon = const_cast<FieldType&>(dimensionedInternalField());
+
+    // only set the values if the weights are < 1 - tolerance
+    forAll(weights, faceI)
+    {
+        scalar w = weights[faceI];
+
+        if (w < 1.0 - 1e-6)
+        {
+            label cellI = patch().faceCells()[faceI];
+
+            G[cellI] = w*G[cellI] + (1.0 - w)*G0[cellI];
+            epsilon[cellI] = w*epsilon[cellI] + (1.0 - w)*epsilon0[cellI];
+        }
+    }
+
+    fvPatchField<scalar>::updateCoeffs();
+}
+
+
+void epsilonWallFunctionFvPatchScalarField::manipulateMatrix
+(
+    fvMatrix<scalar>& matrix
+)
+{
+    if (manipulatedMatrix())
+    {
+        return;
+    }
+
+    matrix.setValues(patch().faceCells(), patchInternalField());
+
+    fvPatchField<scalar>::manipulateMatrix(matrix);
+}
+
+
+void epsilonWallFunctionFvPatchScalarField::manipulateMatrix
+(
+    fvMatrix<scalar>& matrix,
+    const Field<scalar>& weights
+)
+{
+    if (manipulatedMatrix())
+    {
+        return;
+    }
+
+    // filter weights so that we only apply the constraint where the
+    // weight > SMALL
+    DynamicList<label> constraintCells(weights.size());
+    DynamicList<scalar> constraintEpsilon(weights.size());
+    const labelUList& faceCells = patch().faceCells();
+
+    const DimensionedField<scalar, volMesh>& epsilon
+        = dimensionedInternalField();
+
+    label nConstrainedCells = 0;
+
+
+    forAll(weights, faceI)
+    {
+        // only set the values if the weights are < 1 - tolerance
+        if (weights[faceI] < (1.0 - 1e-6))
+        {
+            nConstrainedCells++;
+
+            label cellI = faceCells[faceI];
+
+            constraintCells.append(cellI);
+            constraintEpsilon.append(epsilon[cellI]);
+        }
+    }
+
+    if (debug)
+    {
+        Pout<< "Patch: " << patch().name()
+            << ": number of constrained cells = " << nConstrainedCells
+            << " out of " << patch().size()
+            << endl;
+    }
+
+    matrix.setValues
+    (
+        constraintCells,
+        scalarField(constraintEpsilon.xfer())
+    );
+
+    fvPatchField<scalar>::manipulateMatrix(matrix);
 }
 
 
 void epsilonWallFunctionFvPatchScalarField::write(Ostream& os) const
 {
-    fixedInternalValueFvPatchField<scalar>::write(os);
+    fixedValueFvPatchField<scalar>::write(os);
     writeLocalEntries(os);
     writeEntry("value", os);
 }
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.H b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.H
index 9ff7b7acb5..61de72a385 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.H
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -71,7 +71,7 @@ SourceFiles
 #ifndef epsilonWallFunctionFvPatchScalarField_H
 #define epsilonWallFunctionFvPatchScalarField_H
 
-#include "fixedInternalValueFvPatchField.H"
+#include "fixedValueFvPatchField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -80,13 +80,15 @@ namespace Foam
 namespace incompressible
 {
 
+class turbulenceModel;
+
 /*---------------------------------------------------------------------------*\
            Class epsilonWallFunctionFvPatchScalarField Declaration
 \*---------------------------------------------------------------------------*/
 
 class epsilonWallFunctionFvPatchScalarField
 :
-    public fixedInternalValueFvPatchField<scalar>
+    public fixedValueFvPatchField<scalar>
 {
 protected:
 
@@ -101,6 +103,21 @@ protected:
         //- E coefficient
         scalar E_;
 
+        //- Local copy of turbulence G field
+        scalarField G_;
+
+        //- Local copy of turbulence epsilon field
+        scalarField epsilon_;
+
+        //- Initialised flag
+        bool initialised_;
+
+        //- Master patch ID
+        label master_;
+
+        //- List of averaging corner weights
+        List<List<scalar> > cornerWeights_;
+
 
     // Protected Member Functions
 
@@ -110,6 +127,44 @@ protected:
         //- Write local wall function variables
         virtual void writeLocalEntries(Ostream&) const;
 
+        //- Set the master patch - master is responsible for updating all
+        //  wall function patches
+        virtual void setMaster();
+
+        //- Create the averaging weights for cells which are bounded by
+        //  multiple wall function faces
+        virtual void createAveragingWeights();
+
+        //- Helper function to return non-const access to an epsilon patch
+        virtual epsilonWallFunctionFvPatchScalarField& epsilonPatch
+        (
+            const label patchI
+        );
+
+        //- Main driver to calculate the turbulence fields
+        virtual void calculateTurbulenceFields
+        (
+            const turbulenceModel& turbulence,
+            scalarField& G0,
+            scalarField& epsilon0
+        );
+
+        //- Calculate the epsilon and G
+        virtual void calculate
+        (
+            const turbulenceModel& turbulence,
+            const List<scalar>& cornerWeights,
+            const fvPatch& patch,
+            scalarField& G,
+            scalarField& epsilon
+        );
+
+        //- Return non-const access to the master patch ID
+        virtual label& master()
+        {
+            return master_;
+        }
+
 
 public:
 
@@ -179,16 +234,39 @@ public:
             );
         }
 
+    //- Destructor
+    virtual ~epsilonWallFunctionFvPatchScalarField()
+    {}
+
 
     // Member functions
 
+        // Access
+
+            //- Return non-const access to the master's G field
+            scalarField& G(bool init = false);
+
+            //- Return non-const access to the master's epsilon field
+            scalarField& epsilon(bool init = false);
+
+
         // Evaluation functions
 
             //- Update the coefficients associated with the patch field
             virtual void updateCoeffs();
 
-            //- Evaluate the patchField
-            virtual void evaluate(const Pstream::commsTypes);
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs(const scalarField& weights);
+
+            //- Manipulate matrix
+            virtual void manipulateMatrix(fvMatrix<scalar>& matrix);
+
+            //- Manipulate matrix with given weights
+            virtual void manipulateMatrix
+            (
+                fvMatrix<scalar>& matrix,
+                const scalarField& weights
+            );
 
 
         // I-O