diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.C
index bc3bab534d..905a565b0a 100644
--- a/src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.C
@@ -284,6 +284,16 @@ Foam::tmp<Foam::Field<Type> > Foam::cyclicACMIFvPatchField<Type>::snGrad
 }
 
 
+template<class Type>
+void Foam::cyclicACMIFvPatchField<Type>::updateCoeffs()
+{
+    const scalarField& mask = cyclicACMIPatch_.cyclicACMIPatch().mask();
+
+    const fvPatchField<Type>& npf = nonOverlapPatchField();
+    const_cast<fvPatchField<Type>&>(npf).updateCoeffs(mask);
+}
+
+
 template<class Type>
 void Foam::cyclicACMIFvPatchField<Type>::evaluate
 (
@@ -376,6 +386,20 @@ Foam::cyclicACMIFvPatchField<Type>::gradientBoundaryCoeffs() const
 }
 
 
+template<class Type>
+void Foam::cyclicACMIFvPatchField<Type>::manipulateMatrix
+(
+    fvMatrix<Type>& matrix
+)
+{
+    // blend contrubutions from the coupled and non-overlap patches
+    const fvPatchField<Type>& npf = nonOverlapPatchField();
+
+    const scalarField& mask = cyclicACMIPatch_.cyclicACMIPatch().mask();
+    const_cast<fvPatchField<Type>&>(npf).manipulateMatrix(matrix, mask);
+}
+
+
 template<class Type>
 void Foam::cyclicACMIFvPatchField<Type>::write(Ostream& os) const
 {
diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.H
index 9848505be1..91f414d321 100644
--- a/src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.H
+++ b/src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.H
@@ -178,12 +178,34 @@ public:
             //- Return reference to non-overlapping patchField
             const fvPatchField<Type>& nonOverlapPatchField() const;
 
+            //- Update result field based on interface functionality
+            virtual void updateInterfaceMatrix
+            (
+                scalarField& result,
+                const scalarField& psiInternal,
+                const scalarField& coeffs,
+                const direction cmpt,
+                const Pstream::commsTypes commsType
+            ) const;
+
+            //- Update result field based on interface functionality
+            virtual void updateInterfaceMatrix
+            (
+                Field<Type>&,
+                const Field<Type>&,
+                const scalarField&,
+                const Pstream::commsTypes commsType
+            ) const;
+
             //- Return patch-normal gradient
             virtual tmp<Field<Type> > snGrad
             (
                 const scalarField& deltaCoeffs
             ) const;
 
+            //- Update the coefficients associated with the patch field
+            void updateCoeffs();
+
             //- Evaluate the patch field
             virtual void evaluate
             (
@@ -226,24 +248,8 @@ public:
             //  evaluation of the gradient of this patchField
             virtual tmp<Field<Type> > gradientBoundaryCoeffs() const;
 
-            //- Update result field based on interface functionality
-            virtual void updateInterfaceMatrix
-            (
-                scalarField& result,
-                const scalarField& psiInternal,
-                const scalarField& coeffs,
-                const direction cmpt,
-                const Pstream::commsTypes commsType
-            ) const;
-
-            //- Update result field based on interface functionality
-            virtual void updateInterfaceMatrix
-            (
-                Field<Type>&,
-                const Field<Type>&,
-                const scalarField&,
-                const Pstream::commsTypes commsType
-            ) const;
+            //- Manipulate matrix
+            virtual void manipulateMatrix(fvMatrix<Type>& matrix);
 
 
         // Cyclic AMI coupled interface functions
diff --git a/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.C b/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.C
index fa3d6d7a31..7fd83f41b1 100644
--- a/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.C
@@ -42,6 +42,7 @@ Foam::fvPatchField<Type>::fvPatchField
     patch_(p),
     internalField_(iF),
     updated_(false),
+    manipulatedMatrix_(false),
     patchType_(word::null)
 {}
 
@@ -58,6 +59,7 @@ Foam::fvPatchField<Type>::fvPatchField
     patch_(p),
     internalField_(iF),
     updated_(false),
+    manipulatedMatrix_(false),
     patchType_(word::null)
 {}
 
@@ -75,6 +77,7 @@ Foam::fvPatchField<Type>::fvPatchField
     patch_(p),
     internalField_(iF),
     updated_(false),
+    manipulatedMatrix_(false),
     patchType_(ptf.patchType_)
 {}
 
@@ -92,6 +95,7 @@ Foam::fvPatchField<Type>::fvPatchField
     patch_(p),
     internalField_(iF),
     updated_(false),
+    manipulatedMatrix_(false),
     patchType_(dict.lookupOrDefault<word>("patchType", word::null))
 {
     if (dict.found("value"))
@@ -133,6 +137,7 @@ Foam::fvPatchField<Type>::fvPatchField
     patch_(ptf.patch_),
     internalField_(ptf.internalField_),
     updated_(false),
+    manipulatedMatrix_(false),
     patchType_(ptf.patchType_)
 {}
 
@@ -148,6 +153,7 @@ Foam::fvPatchField<Type>::fvPatchField
     patch_(ptf.patch_),
     internalField_(iF),
     updated_(false),
+    manipulatedMatrix_(false),
     patchType_(ptf.patchType_)
 {}
 
@@ -267,6 +273,28 @@ void Foam::fvPatchField<Type>::rmap
 }
 
 
+template<class Type>
+void Foam::fvPatchField<Type>::updateCoeffs()
+{
+    updated_ = true;
+}
+
+
+template<class Type>
+void Foam::fvPatchField<Type>::updateCoeffs(const scalarField& weights)
+{
+    if (!updated_)
+    {
+        updateCoeffs();
+
+        Field<Type>& fld = *this;
+        fld *= weights;
+
+        updated_ = true;
+    }
+}
+
+
 template<class Type>
 void Foam::fvPatchField<Type>::evaluate(const Pstream::commsTypes)
 {
@@ -276,13 +304,25 @@ void Foam::fvPatchField<Type>::evaluate(const Pstream::commsTypes)
     }
 
     updated_ = false;
+    manipulatedMatrix_ = false;
 }
 
 
 template<class Type>
 void Foam::fvPatchField<Type>::manipulateMatrix(fvMatrix<Type>& matrix)
 {
-    // do nothing
+    manipulatedMatrix_ = true;
+}
+
+
+template<class Type>
+void Foam::fvPatchField<Type>::manipulateMatrix
+(
+    fvMatrix<Type>& matrix,
+    const scalarField& weights
+)
+{
+    manipulatedMatrix_ = true;
 }
 
 
diff --git a/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H b/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H
index 6047f1e4af..d3b2ed4e70 100644
--- a/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H
+++ b/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H
@@ -93,6 +93,10 @@ class fvPatchField
         //  the construction of the matrix
         bool updated_;
 
+        //- Update index used so that manipulateMatrix is called only once
+        //  during the construction of the matrix
+        bool manipulatedMatrix_;
+
         //- Optional patch type, used to allow specified boundary conditions
         //  to be applied to constraint patches by providing the constraint
         //  patch type as 'patchType'
@@ -327,6 +331,12 @@ public:
                 return updated_;
             }
 
+            //- Return true if the matrix has already been manipulated
+            bool manipulatedMatrix() const
+            {
+                return manipulatedMatrix_;
+            }
+
 
         // Mapping functions
 
@@ -365,10 +375,12 @@ public:
 
             //- Update the coefficients associated with the patch field
             //  Sets Updated to true
-            virtual void updateCoeffs()
-            {
-                updated_ = true;
-            }
+            virtual void updateCoeffs();
+
+            //- Update the coefficients associated with the patch field
+            //  and apply weight field
+            //  Sets Updated to true
+            virtual void updateCoeffs(const scalarField& weights);
 
             //- Return internal field next to patch as patch field
             virtual tmp<Field<Type> > patchInternalField() const;
@@ -479,6 +491,13 @@ public:
             //- Manipulate matrix
             virtual void manipulateMatrix(fvMatrix<Type>& matrix);
 
+            //- Manipulate matrix with given weights
+            virtual void manipulateMatrix
+            (
+                fvMatrix<Type>& matrix,
+                const scalarField& weights
+            );
+
 
         // I-O
 
diff --git a/src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.C b/src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.C
index 6e249d876f..30ae0adf37 100644
--- a/src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.C
+++ b/src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.C
@@ -342,6 +342,18 @@ Foam::label Foam::cyclicACMIPolyPatch::nonOverlapPatchID() const
                 << exit(FatalError);
         }
 
+        if (nonOverlapPatchID_ < index())
+        {
+            FatalErrorIn("cyclicPolyAMIPatch::neighbPatchID() const")
+                << "Boundary ordering error: " << type()
+                << " patch must be defined prior to its non-overlapping patch"
+                << nl
+                << type() << " patch: " << name() << ", ID:" << index() << nl
+                << "Non-overlap patch: " << nonOverlapPatchName_
+                << ", ID:" << nonOverlapPatchID_ << nl
+                << exit(FatalError);
+        }
+
         const polyPatch& noPp = this->boundaryMesh()[nonOverlapPatchID_];
 
         bool ok = true;
diff --git a/src/turbulenceModels/compressible/RAS/backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.C b/src/turbulenceModels/compressible/RAS/backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.C
index 55434df625..7846bec820 100644
--- a/src/turbulenceModels/compressible/RAS/backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.C
+++ b/src/turbulenceModels/compressible/RAS/backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.C
@@ -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
@@ -25,6 +25,7 @@ License
 
 #include "backwardsCompatibilityWallFunctions.H"
 
+#include "volFields.H"
 #include "calculatedFvPatchField.H"
 #include "alphatWallFunctionFvPatchScalarField.H"
 #include "mutkWallFunctionFvPatchScalarField.H"
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C
index 84762e6beb..d079c811ea 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C
@@ -55,6 +55,67 @@ 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> tmu = turbulence.mu();
+    const scalarField& muw = tmu().boundaryField()[patchI];
+
+    const tmp<volScalarField> tmut = turbulence.mut();
+    const volScalarField& mut = tmut();
+    const scalarField& mutw = mut.boundaryField()[patchI];
+
+    const scalarField& rhow = turbulence.rho().boundaryField()[patchI];
+
+    const fvPatchVectorField& Uw = turbulence.U().boundaryField()[patchI];
+
+    const scalarField magGradUw(mag(Uw.snGrad()));
+
+    // Set epsilon and G
+    forAll(mutw, faceI)
+    {
+        label cellI = patch.faceCells()[faceI];
+
+        scalar yPlus = Cmu25*sqrt(k[cellI])*y[faceI]/muw[faceI]/rhow[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]*muw[faceI]/rhow[faceI]/sqr(y[faceI]);
+        }
+
+        G[cellI] =
+            w
+           *(mutw[faceI] + muw[faceI])
+           *magGradUw[faceI]
+           *Cmu25*sqrt(k[cellI])
+           /(kappa_*y[faceI]);
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 epsilonLowReWallFunctionFvPatchScalarField::
@@ -119,84 +180,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> tmu = turbulence.mu();
-    const scalarField& muw = tmu().boundaryField()[patchI];
-
-    const tmp<volScalarField> tmut = turbulence.mut();
-    const volScalarField& mut = tmut();
-    const scalarField& mutw = mut.boundaryField()[patchI];
-
-    const scalarField& rhow = turbulence.rho().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(mutw, faceI)
-    {
-        label faceCellI = patch().faceCells()[faceI];
-
-        scalar yPlus = Cmu25*sqrt(k[faceCellI])*y[faceI]/muw[faceI]/rhow[faceI];
-
-        if (yPlus > yPlusLam_)
-        {
-            epsilon[faceCellI] = Cmu75*pow(k[faceCellI], 1.5)/(kappa_*y[faceI]);
-        }
-        else
-        {
-            epsilon[faceCellI] =
-                2.0*k[faceCellI]*muw[faceI]/rhow[faceI]/sqr(y[faceI]);
-        }
-
-        G[faceCellI] =
-            (mutw[faceI] + muw[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/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.H b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.H
index 61434c29b5..c2f7cf178f 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.H
+++ b/src/turbulenceModels/compressible/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:
 
@@ -165,14 +175,6 @@ public:
                 new epsilonLowReWallFunctionFvPatchScalarField(*this, iF)
             );
         }
-
-
-    // Member functions
-
-        // Evaluation functions
-
-            //- Update the coefficients associated with the patch field
-            virtual void updateCoeffs();
 };
 
 
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C
index dc6296f2ed..5e78e48c97 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C
@@ -26,10 +26,10 @@ License
 #include "epsilonWallFunctionFvPatchScalarField.H"
 #include "compressible/turbulenceModel/turbulenceModel.H"
 #include "fvPatchFieldMapper.H"
+#include "fvMatrix.H"
 #include "volFields.H"
-#include "addToRunTimeSelectionTable.H"
-#include "mutWallFunctionFvPatchScalarField.H"
 #include "wallFvPatch.H"
+#include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -62,6 +62,193 @@ void epsilonWallFunctionFvPatchScalarField::writeLocalEntries(Ostream& os) const
 }
 
 
+void epsilonWallFunctionFvPatchScalarField::setMaster()
+{
+    if (master_ != -1)
+    {
+        return;
+    }
+
+    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 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> tmu = turbulence.mu();
+    const scalarField& muw = tmu().boundaryField()[patchI];
+
+    const tmp<volScalarField> tmut = turbulence.mut();
+    const volScalarField& mut = tmut();
+    const scalarField& mutw = mut.boundaryField()[patchI];
+
+    const fvPatchVectorField& Uw = turbulence.U().boundaryField()[patchI];
+
+    const scalarField magGradUw(mag(Uw.snGrad()));
+
+    // Set epsilon and G
+    forAll(mutw, faceI)
+    {
+        label cellI = patch.faceCells()[faceI];
+
+        scalar w = cornerWeights[faceI];
+
+        epsilon[cellI] = w*Cmu75*pow(k[cellI], 1.5)/(kappa_*y[faceI]);
+
+        G[cellI] =
+            w
+           *(mutw[faceI] + muw[faceI])
+           *magGradUw[faceI]
+           *Cmu25*sqrt(k[cellI])
+           /(kappa_*y[faceI]);
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
@@ -70,10 +257,15 @@ epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    fixedInternalValueFvPatchField<scalar>(p, iF),
+    fixedValueFvPatchField<scalar>(p, iF),
     Cmu_(0.09),
     kappa_(0.41),
-    E_(9.8)
+    E_(9.8),
+    G_(),
+    epsilon_(),
+    initialised_(false),
+    master_(-1),
+    cornerWeights_()
 {
     checkType();
 }
@@ -87,10 +279,15 @@ epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
     const fvPatchFieldMapper& mapper
 )
 :
-    fixedInternalValueFvPatchField<scalar>(ptf, p, iF, mapper),
+    fixedValueFvPatchField<scalar>(ptf, p, iF, mapper),
     Cmu_(ptf.Cmu_),
     kappa_(ptf.kappa_),
-    E_(ptf.E_)
+    E_(ptf.E_),
+    G_(),
+    epsilon_(),
+    initialised_(false),
+    master_(-1),
+    cornerWeights_()
 {
     checkType();
 }
@@ -103,10 +300,15 @@ epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
     const dictionary& dict
 )
 :
-    fixedInternalValueFvPatchField<scalar>(p, iF, 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))
+    E_(dict.lookupOrDefault<scalar>("E", 9.8)),
+    G_(),
+    epsilon_(),
+    initialised_(false),
+    master_(-1),
+    cornerWeights_()
 {
     checkType();
 }
@@ -117,10 +319,15 @@ epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
     const epsilonWallFunctionFvPatchScalarField& ewfpsf
 )
 :
-    fixedInternalValueFvPatchField<scalar>(ewfpsf),
+    fixedValueFvPatchField<scalar>(ewfpsf),
     Cmu_(ewfpsf.Cmu_),
     kappa_(ewfpsf.kappa_),
-    E_(ewfpsf.E_)
+    E_(ewfpsf.E_),
+    G_(),
+    epsilon_(),
+    initialised_(false),
+    master_(-1),
+    cornerWeights_()
 {
     checkType();
 }
@@ -132,10 +339,15 @@ epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    fixedInternalValueFvPatchField<scalar>(ewfpsf, iF),
+    fixedValueFvPatchField<scalar>(ewfpsf, iF),
     Cmu_(ewfpsf.Cmu_),
     kappa_(ewfpsf.kappa_),
-    E_(ewfpsf.E_)
+    E_(ewfpsf.E_),
+    G_(),
+    epsilon_(),
+    initialised_(false),
+    master_(-1),
+    cornerWeights_()
 {
     checkType();
 }
@@ -143,6 +355,38 @@ epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
 
 // * * * * * * * * * * * * * * * 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())
@@ -150,76 +394,168 @@ void epsilonWallFunctionFvPatchScalarField::updateCoeffs()
         return;
     }
 
-    const label patchI = patch().index();
-
     const turbulenceModel& turbulence =
-        db().lookupObject<turbulenceModel>("turbulenceModel");
+        db().lookupObject<turbulenceModel>(turbulenceModel::typeName);
 
-    const scalar Cmu25 = pow025(Cmu_);
-    const scalar Cmu75 = pow(Cmu_, 0.75);
+    setMaster();
 
-    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 scalarField& muw = turbulence.mu().boundaryField()[patchI];
-
-    const tmp<volScalarField> tmut = turbulence.mut();
-    const volScalarField& mut = tmut();
-    const scalarField& mutw = mut.boundaryField()[patchI];
-
-    const fvPatchVectorField& Uw = turbulence.U().boundaryField()[patchI];
-
-    const scalarField magGradUw(mag(Uw.snGrad()));
-
-    // Set epsilon and G
-    forAll(mutw, faceI)
+    if (patch().index() == master_)
     {
-        label faceCellI = patch().faceCells()[faceI];
-
-        epsilon[faceCellI] = Cmu75*pow(k[faceCellI], 1.5)/(kappa_*y[faceI]);
-
-        G[faceCellI] =
-            (mutw[faceI] + muw[faceI])
-           *magGradUw[faceI]
-           *Cmu25*sqrt(k[faceCellI])
-           /(kappa_*y[faceI]);
+        createAveragingWeights();
+        calculateTurbulenceFields(turbulence, G(true), epsilon(true));
     }
 
-    fixedInternalValueFvPatchField<scalar>::updateCoeffs();
+    const scalarField& G0 = this->G();
+    const scalarField& epsilon0 = this->epsilon();
 
-    // TODO: perform averaging for cells sharing more than one boundary face
+    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::evaluate
+void epsilonWallFunctionFvPatchScalarField::updateCoeffs
 (
-    const Pstream::commsTypes commsType
+    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/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.H b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.H
index 9975348d0b..9eefc78620 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.H
+++ b/src/turbulenceModels/compressible/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 compressibleEpsilonWallFunctionFvPatchScalarField_H
 #define compressibleEpsilonWallFunctionFvPatchScalarField_H
 
-#include "fixedInternalValueFvPatchField.H"
+#include "fixedValueFvPatchField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -80,13 +80,15 @@ namespace Foam
 namespace compressible
 {
 
+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:
 
@@ -180,15 +235,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
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
index a826fdf5e3..fff165e15e 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
@@ -26,10 +26,11 @@ License
 #include "omegaWallFunctionFvPatchScalarField.H"
 #include "compressible/turbulenceModel/turbulenceModel.H"
 #include "fvPatchFieldMapper.H"
+#include "fvMatrix.H"
 #include "volFields.H"
-#include "addToRunTimeSelectionTable.H"
-#include "mutWallFunctionFvPatchScalarField.H"
 #include "wallFvPatch.H"
+#include "mutWallFunctionFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -63,6 +64,198 @@ void omegaWallFunctionFvPatchScalarField::writeLocalEntries(Ostream& os) const
 }
 
 
+void omegaWallFunctionFvPatchScalarField::setMaster()
+{
+    if (master_ != -1)
+    {
+        return;
+    }
+
+    const volScalarField& omega =
+        static_cast<const volScalarField&>(this->dimensionedInternalField());
+
+    const volScalarField::GeometricBoundaryField& bf = omega.boundaryField();
+
+    label master = -1;
+    forAll(bf, patchI)
+    {
+        if (isA<omegaWallFunctionFvPatchScalarField>(bf[patchI]))
+        {
+            omegaWallFunctionFvPatchScalarField& epf = omegaPatch(patchI);
+
+            if (master == -1)
+            {
+                master = patchI;
+            }
+
+            epf.master() = master;
+        }
+    }
+}
+
+
+void omegaWallFunctionFvPatchScalarField::createAveragingWeights()
+{
+    if (initialised_)
+    {
+        return;
+    }
+
+    const volScalarField& omega =
+        static_cast<const volScalarField&>(this->dimensionedInternalField());
+
+    const volScalarField::GeometricBoundaryField& bf = omega.boundaryField();
+
+    const fvMesh& mesh = omega.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> omegaPatches(bf.size());
+    forAll(bf, patchI)
+    {
+        if (isA<omegaWallFunctionFvPatchScalarField>(bf[patchI]))
+        {
+            omegaPatches.append(patchI);
+
+            const labelUList& faceCells = bf[patchI].patch().faceCells();
+            forAll(faceCells, i)
+            {
+                label cellI = faceCells[i];
+                weights[cellI]++;
+            }
+        }
+    }
+
+    cornerWeights_.setSize(bf.size());
+    forAll(omegaPatches, i)
+    {
+        label patchI = omegaPatches[i];
+        const fvPatchField& wf = weights.boundaryField()[patchI];
+        cornerWeights_[patchI] = 1.0/wf.patchInternalField();
+    }
+
+    G_.setSize(dimensionedInternalField().size(), 0.0);
+    omega_.setSize(dimensionedInternalField().size(), 0.0);
+
+    initialised_ = true;
+}
+
+
+omegaWallFunctionFvPatchScalarField&
+omegaWallFunctionFvPatchScalarField::omegaPatch(const label patchI)
+{
+    const volScalarField& omega =
+        static_cast<const volScalarField&>(this->dimensionedInternalField());
+
+    const volScalarField::GeometricBoundaryField& bf = omega.boundaryField();
+
+    const omegaWallFunctionFvPatchScalarField& epf =
+        refCast<const omegaWallFunctionFvPatchScalarField>(bf[patchI]);
+
+    return const_cast<omegaWallFunctionFvPatchScalarField&>(epf);
+}
+
+
+void omegaWallFunctionFvPatchScalarField::calculateTurbulenceFields
+(
+    const turbulenceModel& turbulence,
+    scalarField& G0,
+    scalarField& omega0
+)
+{
+    // accumulate all of the G and omega contributions
+    forAll(cornerWeights_, patchI)
+    {
+        if (!cornerWeights_[patchI].empty())
+        {
+            omegaWallFunctionFvPatchScalarField& epf = omegaPatch(patchI);
+
+            const List<scalar>& w = cornerWeights_[patchI];
+
+            epf.calculate(turbulence, w, epf.patch(), G0, omega0);
+        }
+    }
+
+    // apply zero-gradient condition for omega
+    forAll(cornerWeights_, patchI)
+    {
+        if (!cornerWeights_[patchI].empty())
+        {
+            omegaWallFunctionFvPatchScalarField& epf = omegaPatch(patchI);
+
+            epf == scalarField(omega0, epf.patch().faceCells());
+        }
+    }
+}
+
+
+void omegaWallFunctionFvPatchScalarField::calculate
+(
+    const turbulenceModel& turbulence,
+    const List<scalar>& cornerWeights,
+    const fvPatch& patch,
+    scalarField& G,
+    scalarField& omega
+)
+{
+    const label patchI = patch.index();
+
+    const scalarField& y = turbulence.y()[patchI];
+
+    const scalar Cmu25 = pow025(Cmu_);
+
+    const tmp<volScalarField> tk = turbulence.k();
+    const volScalarField& k = tk();
+
+    const scalarField& rhow = turbulence.rho().boundaryField()[patchI];
+
+    const tmp<volScalarField> tmu = turbulence.mu();
+    const scalarField& muw = tmu().boundaryField()[patchI];
+
+    const tmp<volScalarField> tmut = turbulence.mut();
+    const volScalarField& mut = tmut();
+    const scalarField& mutw = mut.boundaryField()[patchI];
+
+    const fvPatchVectorField& Uw = turbulence.U().boundaryField()[patchI];
+
+    const scalarField magGradUw(mag(Uw.snGrad()));
+
+    // Set omega and G
+    forAll(mutw, faceI)
+    {
+        label cellI = patch.faceCells()[faceI];
+
+        scalar w = cornerWeights[faceI];
+
+        scalar omegaVis = 6.0*muw[faceI]/(rhow[faceI]*beta1_*sqr(y[faceI]));
+
+        scalar omegaLog = sqrt(k[cellI])/(Cmu25*kappa_*y[faceI]);
+
+        omega[cellI] = w*sqrt(sqr(omegaVis) + sqr(omegaLog));
+
+        G[cellI] =
+            w
+           *(mutw[faceI] + muw[faceI])
+           *magGradUw[faceI]
+           *Cmu25*sqrt(k[cellI])
+           /(kappa_*y[faceI]);
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
@@ -71,13 +264,17 @@ omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    fixedInternalValueFvPatchField<scalar>(p, iF),
+    fixedValueFvPatchField<scalar>(p, iF),
     Cmu_(0.09),
     kappa_(0.41),
     E_(9.8),
     beta1_(0.075),
-    yPlusLam_(mutWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_))
-
+    yPlusLam_(mutWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_)),
+    G_(),
+    omega_(),
+    initialised_(false),
+    master_(-1),
+    cornerWeights_()
 {
     checkType();
 }
@@ -91,12 +288,17 @@ omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
     const fvPatchFieldMapper& mapper
 )
 :
-    fixedInternalValueFvPatchField<scalar>(ptf, p, iF, mapper),
+    fixedValueFvPatchField<scalar>(ptf, p, iF, mapper),
     Cmu_(ptf.Cmu_),
     kappa_(ptf.kappa_),
     E_(ptf.E_),
     beta1_(ptf.beta1_),
-    yPlusLam_(ptf.yPlusLam_)
+    yPlusLam_(ptf.yPlusLam_),
+    G_(),
+    omega_(),
+    initialised_(false),
+    master_(-1),
+    cornerWeights_()
 {
     checkType();
 }
@@ -109,12 +311,17 @@ omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
     const dictionary& dict
 )
 :
-    fixedInternalValueFvPatchField<scalar>(p, iF, 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)),
     beta1_(dict.lookupOrDefault<scalar>("beta1", 0.075)),
-    yPlusLam_(mutWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_))
+    yPlusLam_(mutWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_)),
+    G_(),
+    omega_(),
+    initialised_(false),
+    master_(-1),
+    cornerWeights_()
 {
     checkType();
 }
@@ -125,12 +332,17 @@ omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
     const omegaWallFunctionFvPatchScalarField& owfpsf
 )
 :
-    fixedInternalValueFvPatchField<scalar>(owfpsf),
+    fixedValueFvPatchField<scalar>(owfpsf),
     Cmu_(owfpsf.Cmu_),
     kappa_(owfpsf.kappa_),
     E_(owfpsf.E_),
     beta1_(owfpsf.beta1_),
-    yPlusLam_(owfpsf.yPlusLam_)
+    yPlusLam_(owfpsf.yPlusLam_),
+    G_(),
+    omega_(),
+    initialised_(false),
+    master_(-1),
+    cornerWeights_()
 {
     checkType();
 }
@@ -142,13 +354,17 @@ omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    fixedInternalValueFvPatchField<scalar>(owfpsf, iF),
+    fixedValueFvPatchField<scalar>(owfpsf, iF),
     Cmu_(owfpsf.Cmu_),
     kappa_(owfpsf.kappa_),
     E_(owfpsf.E_),
     beta1_(owfpsf.beta1_),
-    yPlusLam_(owfpsf.yPlusLam_)
-
+    yPlusLam_(owfpsf.yPlusLam_),
+    G_(),
+    omega_(),
+    initialised_(false),
+    master_(-1),
+    cornerWeights_()
 {
     checkType();
 }
@@ -156,6 +372,38 @@ omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+scalarField& omegaWallFunctionFvPatchScalarField::G(bool init)
+{
+    if (patch().index() == master_)
+    {
+        if (init)
+        {
+            G_ = 0.0;
+        }
+
+        return G_;
+    }
+
+    return omegaPatch(master_).G();
+}
+
+
+scalarField& omegaWallFunctionFvPatchScalarField::omega(bool init)
+{
+    if (patch().index() == master_)
+    {
+        if (init)
+        {
+            omega_ = 0.0;
+        }
+
+        return omega_;
+    }
+
+    return omegaPatch(master_).omega(init);
+}
+
+
 void omegaWallFunctionFvPatchScalarField::updateCoeffs()
 {
     if (updated())
@@ -163,68 +411,168 @@ void omegaWallFunctionFvPatchScalarField::updateCoeffs()
         return;
     }
 
-    const label patchI = patch().index();
-
     const turbulenceModel& turbulence =
-        db().lookupObject<turbulenceModel>("turbulenceModel");
-    const scalarField& y = turbulence.y()[patch().index()];
+        db().lookupObject<turbulenceModel>(turbulenceModel::typeName);
 
-    const scalar Cmu25 = pow025(Cmu_);
+    setMaster();
 
-    volScalarField& G = const_cast<volScalarField&>
-        (
-            db().lookupObject<volScalarField>
-            (
-                turbulence.GName()
-            )
-        );
-
-    DimensionedField<scalar, volMesh>& omega =
-        const_cast<DimensionedField<scalar, volMesh>&>
-        (
-            dimensionedInternalField()
-        );
-
-    const tmp<volScalarField> tk = turbulence.k();
-    const volScalarField& k = tk();
-
-    const scalarField& rhow = turbulence.rho().boundaryField()[patchI];
-
-    const scalarField& muw = turbulence.mu().boundaryField()[patchI];
-
-    const tmp<volScalarField> tmut = turbulence.mut();
-    const volScalarField& mut = tmut();
-    const scalarField& mutw = mut.boundaryField()[patchI];
-
-    const fvPatchVectorField& Uw = turbulence.U().boundaryField()[patchI];
-
-    const scalarField magGradUw(mag(Uw.snGrad()));
-
-    // Set omega and G
-    forAll(mutw, faceI)
+    if (patch().index() == master_)
     {
-        label faceCellI = patch().faceCells()[faceI];
-
-        scalar omegaVis = 6.0*muw[faceI]/(rhow[faceI]*beta1_*sqr(y[faceI]));
-        scalar omegaLog = sqrt(k[faceCellI])/(Cmu25*kappa_*y[faceI]);
-        omega[faceCellI] = sqrt(sqr(omegaVis) + sqr(omegaLog));
-
-        G[faceCellI] =
-            (mutw[faceI] + muw[faceI])
-           *magGradUw[faceI]
-           *Cmu25*sqrt(k[faceCellI])
-           /(kappa_*y[faceI]);
+        createAveragingWeights();
+        calculateTurbulenceFields(turbulence, G(true), omega(true));
     }
 
-    fixedInternalValueFvPatchField<scalar>::updateCoeffs();
+    const scalarField& G0 = this->G();
+    const scalarField& omega0 = this->omega();
 
-    // TODO: perform averaging for cells sharing more than one boundary face
+    typedef DimensionedField<scalar, volMesh> FieldType;
+
+    FieldType& G =
+        const_cast<FieldType&>
+        (
+            db().lookupObject<FieldType>(turbulence.GName())
+        );
+
+    FieldType& omega = const_cast<FieldType&>(dimensionedInternalField());
+
+    forAll(*this, faceI)
+    {
+        label cellI = patch().faceCells()[faceI];
+
+        G[cellI] = G0[cellI];
+        omega[cellI] = omega0[cellI];
+    }
+
+    fvPatchField<scalar>::updateCoeffs();
+}
+
+
+void omegaWallFunctionFvPatchScalarField::updateCoeffs
+(
+    const scalarField& weights
+)
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const turbulenceModel& turbulence =
+        db().lookupObject<turbulenceModel>(turbulenceModel::typeName);
+
+    setMaster();
+
+    if (patch().index() == master_)
+    {
+        createAveragingWeights();
+        calculateTurbulenceFields(turbulence, G(true), omega(true));
+    }
+
+    const scalarField& G0 = this->G();
+    const scalarField& omega0 = this->omega();
+
+    typedef DimensionedField<scalar, volMesh> FieldType;
+
+    FieldType& G =
+        const_cast<FieldType&>
+        (
+            db().lookupObject<FieldType>(turbulence.GName())
+        );
+
+    FieldType& omega = 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];
+            omega[cellI] = w*omega[cellI] + (1.0 - w)*omega0[cellI];
+        }
+    }
+
+    fvPatchField<scalar>::updateCoeffs();
+}
+
+
+void omegaWallFunctionFvPatchScalarField::manipulateMatrix
+(
+    fvMatrix<scalar>& matrix
+)
+{
+    if (manipulatedMatrix())
+    {
+        return;
+    }
+
+    matrix.setValues(patch().faceCells(), patchInternalField());
+
+    fvPatchField<scalar>::manipulateMatrix(matrix);
+}
+
+
+void omegaWallFunctionFvPatchScalarField::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> constraintomega(weights.size());
+    const labelUList& faceCells = patch().faceCells();
+
+    const DimensionedField<scalar, volMesh>& omega
+        = 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);
+            constraintomega.append(omega[cellI]);
+        }
+    }
+
+    if (debug)
+    {
+        Pout<< "Patch: " << patch().name()
+            << ": number of constrained cells = " << nConstrainedCells
+            << " out of " << patch().size()
+            << endl;
+    }
+
+    matrix.setValues
+    (
+        constraintCells,
+        scalarField(constraintomega.xfer())
+    );
+
+    fvPatchField<scalar>::manipulateMatrix(matrix);
 }
 
 
 void omegaWallFunctionFvPatchScalarField::write(Ostream& os) const
 {
-    fixedInternalValueFvPatchField<scalar>::write(os);
+    fixedValueFvPatchField<scalar>::write(os);
     writeLocalEntries(os);
     writeEntry("value", os);
 }
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.H b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.H
index 69fa8ee9cd..9db8e32a8d 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.H
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.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
@@ -76,7 +76,7 @@ SourceFiles
 #ifndef compressibleOmegaWallFunctionFvPatchScalarField_H
 #define compressibleOmegaWallFunctionFvPatchScalarField_H
 
-#include "fixedInternalValueFvPatchField.H"
+#include "fixedValueFvPatchField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -85,13 +85,15 @@ namespace Foam
 namespace compressible
 {
 
+class turbulenceModel;
+
 /*---------------------------------------------------------------------------*\
-           Class omegaWallFunctionFvPatchScalarField Declaration
+             Class omegaWallFunctionFvPatchScalarField Declaration
 \*---------------------------------------------------------------------------*/
 
 class omegaWallFunctionFvPatchScalarField
 :
-    public fixedInternalValueFvPatchField<scalar>
+    public fixedValueFvPatchField<scalar>
 {
 protected:
 
@@ -112,6 +114,21 @@ protected:
         //- Y+ at the edge of the laminar sublayer
         scalar yPlusLam_;
 
+        //- Local copy of turbulence G field
+        scalarField G_;
+
+        //- Local copy of turbulence omega field
+        scalarField omega_;
+
+        //- Initialised flag
+        bool initialised_;
+
+        //- Master patch ID
+        label master_;
+
+        //- List of averaging corner weights
+        List<List<scalar> > cornerWeights_;
+
 
     // Protected Member Functions
 
@@ -121,6 +138,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 omega patch
+        virtual omegaWallFunctionFvPatchScalarField& omegaPatch
+        (
+            const label patchI
+        );
+
+        //- Main driver to calculate the turbulence fields
+        virtual void calculateTurbulenceFields
+        (
+            const turbulenceModel& turbulence,
+            scalarField& G0,
+            scalarField& omega0
+        );
+
+        //- Calculate the omega and G
+        virtual void calculate
+        (
+            const turbulenceModel& turbulence,
+            const List<scalar>& cornerWeights,
+            const fvPatch& patch,
+            scalarField& G,
+            scalarField& omega
+        );
+
+        //- Return non-const access to the master patch ID
+        virtual label& master()
+        {
+            return master_;
+        }
+
 
 public:
 
@@ -193,11 +248,33 @@ public:
 
     // 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 omega field
+            scalarField& omega(bool init = false);
+
+
         // Evaluation functions
 
             //- Update the coefficients associated with the patch field
             virtual void updateCoeffs();
 
+            //- 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
 
diff --git a/src/turbulenceModels/incompressible/RAS/backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.C b/src/turbulenceModels/incompressible/RAS/backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.C
index 98514bb5be..df20a315f8 100644
--- a/src/turbulenceModels/incompressible/RAS/backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.C
+++ b/src/turbulenceModels/incompressible/RAS/backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.C
@@ -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
@@ -25,6 +25,7 @@ License
 
 #include "backwardsCompatibilityWallFunctions.H"
 
+#include "volFields.H"
 #include "calculatedFvPatchField.H"
 #include "nutkWallFunctionFvPatchScalarField.H"
 #include "nutLowReWallFunctionFvPatchScalarField.H"
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
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
index a78e0bdb84..554c18542d 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
@@ -26,11 +26,11 @@ License
 #include "omegaWallFunctionFvPatchScalarField.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"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -64,127 +64,159 @@ void omegaWallFunctionFvPatchScalarField::writeLocalEntries(Ostream& os) const
 }
 
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedInternalValueFvPatchField<scalar>(p, iF),
-    Cmu_(0.09),
-    kappa_(0.41),
-    E_(9.8),
-    beta1_(0.075),
-    yPlusLam_(nutkWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_))
+void omegaWallFunctionFvPatchScalarField::setMaster()
 {
-    checkType();
-}
-
-
-omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
-(
-    const omegaWallFunctionFvPatchScalarField& 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_),
-    beta1_(ptf.beta1_),
-    yPlusLam_(ptf.yPlusLam_)
-{
-    checkType();
-}
-
-
-omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
-(
-    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)),
-    beta1_(dict.lookupOrDefault<scalar>("beta1", 0.075)),
-    yPlusLam_(nutkWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_))
-{
-    checkType();
-}
-
-
-omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
-(
-    const omegaWallFunctionFvPatchScalarField& owfpsf
-)
-:
-    fixedInternalValueFvPatchField<scalar>(owfpsf),
-    Cmu_(owfpsf.Cmu_),
-    kappa_(owfpsf.kappa_),
-    E_(owfpsf.E_),
-    beta1_(owfpsf.beta1_),
-    yPlusLam_(owfpsf.yPlusLam_)
-{
-    checkType();
-}
-
-
-omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
-(
-    const omegaWallFunctionFvPatchScalarField& owfpsf,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedInternalValueFvPatchField<scalar>(owfpsf, iF),
-    Cmu_(owfpsf.Cmu_),
-    kappa_(owfpsf.kappa_),
-    E_(owfpsf.E_),
-    beta1_(owfpsf.beta1_),
-    yPlusLam_(owfpsf.yPlusLam_)
-{
-    checkType();
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void omegaWallFunctionFvPatchScalarField::updateCoeffs()
-{
-    if (updated())
+    if (master_ != -1)
     {
         return;
     }
 
-    const label patchI = patch().index();
+    const volScalarField& omega =
+        static_cast<const volScalarField&>(this->dimensionedInternalField());
+
+    const volScalarField::GeometricBoundaryField& bf = omega.boundaryField();
+
+    label master = -1;
+    forAll(bf, patchI)
+    {
+        if (isA<omegaWallFunctionFvPatchScalarField>(bf[patchI]))
+        {
+            omegaWallFunctionFvPatchScalarField& epf = omegaPatch(patchI);
+
+            if (master == -1)
+            {
+                master = patchI;
+            }
+
+            epf.master() = master;
+        }
+    }
+}
+
+
+void omegaWallFunctionFvPatchScalarField::createAveragingWeights()
+{
+    if (initialised_)
+    {
+        return;
+    }
+
+    const volScalarField& omega =
+        static_cast<const volScalarField&>(this->dimensionedInternalField());
+
+    const volScalarField::GeometricBoundaryField& bf = omega.boundaryField();
+
+    const fvMesh& mesh = omega.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> omegaPatches(bf.size());
+    forAll(bf, patchI)
+    {
+        if (isA<omegaWallFunctionFvPatchScalarField>(bf[patchI]))
+        {
+            omegaPatches.append(patchI);
+
+            const labelUList& faceCells = bf[patchI].patch().faceCells();
+            forAll(faceCells, i)
+            {
+                label cellI = faceCells[i];
+                weights[cellI]++;
+            }
+        }
+    }
+
+    cornerWeights_.setSize(bf.size());
+    forAll(omegaPatches, i)
+    {
+        label patchI = omegaPatches[i];
+        const fvPatchField& wf = weights.boundaryField()[patchI];
+        cornerWeights_[patchI] = 1.0/wf.patchInternalField();
+    }
+
+    G_.setSize(dimensionedInternalField().size(), 0.0);
+    omega_.setSize(dimensionedInternalField().size(), 0.0);
+
+    initialised_ = true;
+}
+
+
+omegaWallFunctionFvPatchScalarField&
+omegaWallFunctionFvPatchScalarField::omegaPatch(const label patchI)
+{
+    const volScalarField& omega =
+        static_cast<const volScalarField&>(this->dimensionedInternalField());
+
+    const volScalarField::GeometricBoundaryField& bf = omega.boundaryField();
+
+    const omegaWallFunctionFvPatchScalarField& epf =
+        refCast<const omegaWallFunctionFvPatchScalarField>(bf[patchI]);
+
+    return const_cast<omegaWallFunctionFvPatchScalarField&>(epf);
+}
+
+
+void omegaWallFunctionFvPatchScalarField::calculateTurbulenceFields
+(
+    const turbulenceModel& turbulence,
+    scalarField& G0,
+    scalarField& omega0
+)
+{
+    // accumulate all of the G and omega contributions
+    forAll(cornerWeights_, patchI)
+    {
+        if (!cornerWeights_[patchI].empty())
+        {
+            omegaWallFunctionFvPatchScalarField& epf = omegaPatch(patchI);
+
+            const List<scalar>& w = cornerWeights_[patchI];
+
+            epf.calculate(turbulence, w, epf.patch(), G0, omega0);
+        }
+    }
+
+    // apply zero-gradient condition for omega
+    forAll(cornerWeights_, patchI)
+    {
+        if (!cornerWeights_[patchI].empty())
+        {
+            omegaWallFunctionFvPatchScalarField& epf = omegaPatch(patchI);
+
+            epf == scalarField(omega0, epf.patch().faceCells());
+        }
+    }
+}
+
+
+void omegaWallFunctionFvPatchScalarField::calculate
+(
+    const turbulenceModel& turbulence,
+    const List<scalar>& cornerWeights,
+    const fvPatch& patch,
+    scalarField& G,
+    scalarField& omega
+)
+{
+    const label patchI = patch.index();
 
-    const turbulenceModel& turbulence =
-        db().lookupObject<turbulenceModel>("turbulenceModel");
     const scalarField& y = turbulence.y()[patchI];
 
     const scalar Cmu25 = pow025(Cmu_);
 
-    volScalarField& G =
-        const_cast<volScalarField&>
-        (
-            db().lookupObject<volScalarField>
-            (
-                turbulence.GName()
-            )
-        );
-
-    DimensionedField<scalar, volMesh>& omega =
-        const_cast<DimensionedField<scalar, volMesh>&>
-        (
-            dimensionedInternalField()
-        );
-
     const tmp<volScalarField> tk = turbulence.k();
     const volScalarField& k = tk();
 
@@ -202,30 +234,343 @@ void omegaWallFunctionFvPatchScalarField::updateCoeffs()
     // Set omega and G
     forAll(nutw, faceI)
     {
-        label faceCellI = patch().faceCells()[faceI];
+        label cellI = patch.faceCells()[faceI];
+
+        scalar w = cornerWeights[faceI];
 
         scalar omegaVis = 6.0*nuw[faceI]/(beta1_*sqr(y[faceI]));
 
-        scalar omegaLog = sqrt(k[faceCellI])/(Cmu25*kappa_*y[faceI]);
+        scalar omegaLog = sqrt(k[cellI])/(Cmu25*kappa_*y[faceI]);
 
-        omega[faceCellI] = sqrt(sqr(omegaVis) + sqr(omegaLog));
+        omega[cellI] = w*sqrt(sqr(omegaVis) + sqr(omegaLog));
 
-        G[faceCellI] =
-            (nutw[faceI] + nuw[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
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchField<scalar>(p, iF),
+    Cmu_(0.09),
+    kappa_(0.41),
+    E_(9.8),
+    beta1_(0.075),
+    yPlusLam_(nutkWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_)),
+    G_(),
+    omega_(),
+    initialised_(false),
+    master_(-1),
+    cornerWeights_()
+{
+    checkType();
+}
+
+
+omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
+(
+    const omegaWallFunctionFvPatchScalarField& 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_),
+    beta1_(ptf.beta1_),
+    yPlusLam_(ptf.yPlusLam_),
+    G_(),
+    omega_(),
+    initialised_(false),
+    master_(-1),
+    cornerWeights_()
+{
+    checkType();
+}
+
+
+omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
+(
+    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)),
+    beta1_(dict.lookupOrDefault<scalar>("beta1", 0.075)),
+    yPlusLam_(nutkWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_)),
+    G_(),
+    omega_(),
+    initialised_(false),
+    master_(-1),
+    cornerWeights_()
+{
+    checkType();
+}
+
+
+omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
+(
+    const omegaWallFunctionFvPatchScalarField& owfpsf
+)
+:
+    fixedValueFvPatchField<scalar>(owfpsf),
+    Cmu_(owfpsf.Cmu_),
+    kappa_(owfpsf.kappa_),
+    E_(owfpsf.E_),
+    beta1_(owfpsf.beta1_),
+    yPlusLam_(owfpsf.yPlusLam_),
+    G_(),
+    omega_(),
+    initialised_(false),
+    master_(-1),
+    cornerWeights_()
+{
+    checkType();
+}
+
+
+omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
+(
+    const omegaWallFunctionFvPatchScalarField& owfpsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchField<scalar>(owfpsf, iF),
+    Cmu_(owfpsf.Cmu_),
+    kappa_(owfpsf.kappa_),
+    E_(owfpsf.E_),
+    beta1_(owfpsf.beta1_),
+    yPlusLam_(owfpsf.yPlusLam_),
+    G_(),
+    omega_(),
+    initialised_(false),
+    master_(-1),
+    cornerWeights_()
+{
+    checkType();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+scalarField& omegaWallFunctionFvPatchScalarField::G(bool init)
+{
+    if (patch().index() == master_)
+    {
+        if (init)
+        {
+            G_ = 0.0;
+        }
+
+        return G_;
+    }
+
+    return omegaPatch(master_).G();
+}
+
+
+scalarField& omegaWallFunctionFvPatchScalarField::omega(bool init)
+{
+    if (patch().index() == master_)
+    {
+        if (init)
+        {
+            omega_ = 0.0;
+        }
+
+        return omega_;
+    }
+
+    return omegaPatch(master_).omega(init);
+}
+
+
+void omegaWallFunctionFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const turbulenceModel& turbulence =
+        db().lookupObject<turbulenceModel>(turbulenceModel::typeName);
+
+    setMaster();
+
+    if (patch().index() == master_)
+    {
+        createAveragingWeights();
+        calculateTurbulenceFields(turbulence, G(true), omega(true));
+    }
+
+    const scalarField& G0 = this->G();
+    const scalarField& omega0 = this->omega();
+
+    typedef DimensionedField<scalar, volMesh> FieldType;
+
+    FieldType& G =
+        const_cast<FieldType&>
+        (
+            db().lookupObject<FieldType>(turbulence.GName())
+        );
+
+    FieldType& omega = const_cast<FieldType&>(dimensionedInternalField());
+
+    forAll(*this, faceI)
+    {
+        label cellI = patch().faceCells()[faceI];
+
+        G[cellI] = G0[cellI];
+        omega[cellI] = omega0[cellI];
+    }
+
+    fvPatchField<scalar>::updateCoeffs();
+}
+
+
+void omegaWallFunctionFvPatchScalarField::updateCoeffs
+(
+    const scalarField& weights
+)
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const turbulenceModel& turbulence =
+        db().lookupObject<turbulenceModel>(turbulenceModel::typeName);
+
+    setMaster();
+
+    if (patch().index() == master_)
+    {
+        createAveragingWeights();
+        calculateTurbulenceFields(turbulence, G(true), omega(true));
+    }
+
+    const scalarField& G0 = this->G();
+    const scalarField& omega0 = this->omega();
+
+    typedef DimensionedField<scalar, volMesh> FieldType;
+
+    FieldType& G =
+        const_cast<FieldType&>
+        (
+            db().lookupObject<FieldType>(turbulence.GName())
+        );
+
+    FieldType& omega = 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];
+            omega[cellI] = w*omega[cellI] + (1.0 - w)*omega0[cellI];
+        }
+    }
+
+    fvPatchField<scalar>::updateCoeffs();
+}
+
+
+void omegaWallFunctionFvPatchScalarField::manipulateMatrix
+(
+    fvMatrix<scalar>& matrix
+)
+{
+    if (manipulatedMatrix())
+    {
+        return;
+    }
+
+    matrix.setValues(patch().faceCells(), patchInternalField());
+
+    fvPatchField<scalar>::manipulateMatrix(matrix);
+}
+
+
+void omegaWallFunctionFvPatchScalarField::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> constraintomega(weights.size());
+    const labelUList& faceCells = patch().faceCells();
+
+    const DimensionedField<scalar, volMesh>& omega
+        = 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);
+            constraintomega.append(omega[cellI]);
+        }
+    }
+
+    if (debug)
+    {
+        Pout<< "Patch: " << patch().name()
+            << ": number of constrained cells = " << nConstrainedCells
+            << " out of " << patch().size()
+            << endl;
+    }
+
+    matrix.setValues
+    (
+        constraintCells,
+        scalarField(constraintomega.xfer())
+    );
+
+    fvPatchField<scalar>::manipulateMatrix(matrix);
 }
 
 
 void omegaWallFunctionFvPatchScalarField::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/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.H b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.H
index 9ee99bd6c6..957d0f945e 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.H
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.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
@@ -76,7 +76,7 @@ SourceFiles
 #ifndef omegaWallFunctionFvPatchScalarField_H
 #define omegaWallFunctionFvPatchScalarField_H
 
-#include "fixedInternalValueFvPatchField.H"
+#include "fixedValueFvPatchField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -85,13 +85,15 @@ namespace Foam
 namespace incompressible
 {
 
+class turbulenceModel;
+
 /*---------------------------------------------------------------------------*\
-           Class omegaWallFunctionFvPatchScalarField Declaration
+             Class omegaWallFunctionFvPatchScalarField Declaration
 \*---------------------------------------------------------------------------*/
 
 class omegaWallFunctionFvPatchScalarField
 :
-    public fixedInternalValueFvPatchField<scalar>
+    public fixedValueFvPatchField<scalar>
 {
 protected:
 
@@ -112,6 +114,21 @@ protected:
         //- Y+ at the edge of the laminar sublayer
         scalar yPlusLam_;
 
+        //- Local copy of turbulence G field
+        scalarField G_;
+
+        //- Local copy of turbulence omega field
+        scalarField omega_;
+
+        //- Initialised flag
+        bool initialised_;
+
+        //- Master patch ID
+        label master_;
+
+        //- List of averaging corner weights
+        List<List<scalar> > cornerWeights_;
+
 
     // Protected Member Functions
 
@@ -121,6 +138,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 omega patch
+        virtual omegaWallFunctionFvPatchScalarField& omegaPatch
+        (
+            const label patchI
+        );
+
+        //- Main driver to calculate the turbulence fields
+        virtual void calculateTurbulenceFields
+        (
+            const turbulenceModel& turbulence,
+            scalarField& G0,
+            scalarField& omega0
+        );
+
+        //- Calculate the omega and G
+        virtual void calculate
+        (
+            const turbulenceModel& turbulence,
+            const List<scalar>& cornerWeights,
+            const fvPatch& patch,
+            scalarField& G,
+            scalarField& omega
+        );
+
+        //- Return non-const access to the master patch ID
+        virtual label& master()
+        {
+            return master_;
+        }
+
 
 public:
 
@@ -193,11 +248,33 @@ public:
 
     // 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 omega field
+            scalarField& omega(bool init = false);
+
+
         // Evaluation functions
 
             //- Update the coefficients associated with the patch field
             virtual void updateCoeffs();
 
+            //- 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
 
diff --git a/tutorials/incompressible/pimpleDyMFoam/oscillatingInletACMI2D/system/createBafflesDict b/tutorials/incompressible/pimpleDyMFoam/oscillatingInletACMI2D/system/createBafflesDict
index c3cf5971da..54e2ab822b 100644
--- a/tutorials/incompressible/pimpleDyMFoam/oscillatingInletACMI2D/system/createBafflesDict
+++ b/tutorials/incompressible/pimpleDyMFoam/oscillatingInletACMI2D/system/createBafflesDict
@@ -23,6 +23,9 @@ internalFacesOnly false;
 // Baffles to create.
 baffles
 {
+    // NOTE: cyclicAMI patches MUST BE defined PRIOR to their associted
+    //       blockage patches
+
     ACMI1
     {
         //- Use predefined faceZone to select faces and orientation.
@@ -32,19 +35,6 @@ baffles
         patches
         {
             master
-            {
-                //- Master side patch
-                name            ACMI1_blockage;
-                type            wall;
-            }
-            slave // not used since we're manipulating a boundary patch
-            {
-                //- Slave side patch
-                name            ACMI1_blockage;
-                type            wall;
-            }
-
-            master2
             {
                 //- Master side patch
                 name            ACMI1_couple;
@@ -54,12 +44,26 @@ baffles
                 nonOverlapPatch ACMI1_blockage;
                 transform       noOrdering;
             }
-            slave2 // not used since we're manipulating a boundary patch
+            slave // not used since we're manipulating a boundary patch
             {
                 //- Slave side patch
                 name            ACMI1_couple;
                 type            patch;
             }
+
+            master2
+            {
+                //- Master side patch
+                name            ACMI1_blockage;
+                type            wall;
+            }
+            slave2 // not used since we're manipulating a boundary patch
+            {
+                //- Slave side patch
+                name            ACMI1_blockage;
+                type            wall;
+            }
+
         }
     }
     ACMI2
@@ -71,19 +75,6 @@ baffles
         patches
         {
             master
-            {
-                //- Master side patch
-                name            ACMI2_blockage;
-                type            wall;
-            }
-            slave // not used since we're manipulating a boundary patch
-            {
-                //- Slave side patch
-                name            ACMI2_blockage;
-                type            wall;
-            }
-
-            master2
             {
                 //- Master side patch
                 name            ACMI2_couple;
@@ -93,12 +84,25 @@ baffles
                 nonOverlapPatch ACMI2_blockage;
                 transform       noOrdering;
             }
-            slave2 // not used since we're manipulating a boundary patch
+            slave // not used since we're manipulating a boundary patch
             {
                 //- Slave side patch
                 name            ACMI2_couple;
                 type            patch;
             }
+
+            master2
+            {
+                //- Master side patch
+                name            ACMI2_blockage;
+                type            wall;
+            }
+            slave2 // not used since we're manipulating a boundary patch
+            {
+                //- Slave side patch
+                name            ACMI2_blockage;
+                type            wall;
+            }
         }
     }
 }