Merge branch 'master' into feature/cvMesh

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

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

@@ -65,7 +65,8 @@ autoPtr<refinementSurfaces> createRefinementSurfaces
 (
     const searchableSurfaces& allGeometry,
     const dictionary& surfacesDict,
-    const dictionary& shapeControlDict
+    const dictionary& shapeControlDict,
+    const label gapLevelIncrement
 )
 {
     autoPtr<refinementSurfaces> surfacePtr;
@@ -145,7 +146,7 @@ autoPtr<refinementSurfaces> createRefinementSurfaces
                 globalLevelIncr[surfI] = shapeDict.lookupOrDefault
                 (
                     "gapLevelIncrement",
-                    0
+                    gapLevelIncrement
                 );
             }
             else
@@ -356,7 +357,7 @@ autoPtr<refinementSurfaces> createRefinementSurfaces
                         label levelIncr = regionDict.lookupOrDefault
                         (
                             "gapLevelIncrement",
-                            0
+                            gapLevelIncrement
                         );
                         regionLevelIncr[surfI].insert(regionI, levelIncr);
 
@@ -1031,7 +1032,8 @@ int main(int argc, char *argv[])
             (
                 allGeometry,
                 conformationDict,
-                shapeControlDict
+                shapeControlDict,
+                refineDict.lookupOrDefault("gapLevelIncrement", 0)
             );
     }
     else
@@ -1041,7 +1043,8 @@ int main(int argc, char *argv[])
             new refinementSurfaces
             (
                 allGeometry,
-                refineDict.subDict("refinementSurfaces")
+                refineDict.subDict("refinementSurfaces"),
+                refineDict.lookupOrDefault("gapLevelIncrement", 0)
             )
         );
 
@@ -1052,7 +1055,6 @@ int main(int argc, char *argv[])
     refinementSurfaces& surfaces = surfacesPtr();
 
 
-
     // Checking only?
 
     if (checkGeometry)
@@ -1351,6 +1353,13 @@ int main(int argc, char *argv[])
     const Switch wantSnap(meshDict.lookup("snap"));
     const Switch wantLayers(meshDict.lookup("addLayers"));
 
+    // Refinement parameters
+    const refinementParameters refineParams(refineDict);
+
+    // Snap parameters
+    const snapParameters snapParams(snapDict);
+
+
     if (wantRefine)
     {
         cpuTime timer;
@@ -1364,15 +1373,20 @@ int main(int argc, char *argv[])
             globalToSlavePatch
         );
 
-        // Refinement parameters
-        refinementParameters refineParams(refineDict);
 
         if (!overwrite && !debug)
         {
             const_cast<Time&>(mesh.time())++;
         }
 
-        refineDriver.doRefine(refineDict, refineParams, wantSnap, motionDict);
+        refineDriver.doRefine
+        (
+            refineDict,
+            refineParams,
+            snapParams,
+            wantSnap,
+            motionDict
+        );
 
         writeMesh
         (
@@ -1396,20 +1410,14 @@ int main(int argc, char *argv[])
             globalToSlavePatch
         );
 
-        // Snap parameters
-        snapParameters snapParams(snapDict);
-        // Temporary hack to get access to resolveFeatureAngle
-        scalar curvature;
-        {
-            refinementParameters refineParams(refineDict);
-            curvature = refineParams.curvature();
-        }
-
         if (!overwrite && !debug)
         {
             const_cast<Time&>(mesh.time())++;
         }
 
+        // Use the resolveFeatureAngle from the refinement parameters
+        scalar curvature = refineParams.curvature();
+
         snapDriver.doSnap(snapDict, motionDict, curvature, snapParams);
 
         writeMesh
@@ -1437,17 +1445,12 @@ int main(int argc, char *argv[])
         // Layer addition parameters
         layerParameters layerParams(layerDict, mesh.boundaryMesh());
 
-        //!!! Temporary hack to get access to maxLocalCells
+        // Use the maxLocalCells from the refinement parameters
-        bool preBalance;
+        bool preBalance = returnReduce
-        {
-            refinementParameters refineParams(refineDict);
-
-            preBalance = returnReduce
         (
             (mesh.nCells() >= refineParams.maxLocalCells()),
             orOp<bool>()
         );
-        }
 
 
         if (!overwrite &&  !debug)

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict

@@ -152,6 +152,10 @@ castellatedMeshControls
                 inGroups (meshedPatches);
             }
 
+
+            //- Optional increment (on top of max level) in small gaps
+            //gapLevelIncrement 2;
+
             //- Optional angle to detect small-large cell situation
             //  perpendicular to the surface. Is the angle of face w.r.t.
             //  the local surface normal. Use on flat(ish) surfaces only.
@@ -180,11 +184,17 @@ castellatedMeshControls
     // - used if feature snapping (see snapControls below) is used
     resolveFeatureAngle 30;
 
+    //- Optional increment (on top of max level) in small gaps
+    //gapLevelIncrement 2;
+
+
     // Planar angle:
-    // - used to determine if neighbouring surface normals
+    // - used to determine if surface normals
-    //   are roughly the same so e.g. free-standing baffles can be merged
+    //   are roughly the same or opposite. Used e.g. in proximity refinement
+    //   and to decide when to merge free-standing baffles
+    //
     // If not specified same as resolveFeatureAngle
-    //planarAngle 30;
+    planarAngle 30;
 
 
     // Region-wise refinement
@@ -229,10 +239,8 @@ castellatedMeshControls
     // free-standing zone faces. Not used if there are no faceZones.
     allowFreeStandingZoneFaces true;
 
-
+    //
-    // Optional: whether all baffles get eroded away. WIP. Used for
+    //useTopologicalSnapDetection false;
-    // surface simplification.
-    //allowFreeStandingBaffles false;
 }
 
 // Settings for the snapping.
@@ -462,20 +470,21 @@ meshQualityControls
     // 1 = hex, <= 0 = folded or flattened illegal cell
     minDeterminant 0.001;
 
-    // minFaceWeight (0 -> 0.5)
+    // Relative position of face in relation to cell centres (0.5 for orthogonal
+    // mesh) (0 -> 0.5)
     minFaceWeight 0.05;
 
-    // minVolRatio (0 -> 1)
+    // Volume ratio of neighbouring cells (0 -> 1)
     minVolRatio 0.01;
 
     // must be >0 for Fluent compatibility
     minTriangleTwist -1;
 
-    //- if >0 : preserve single cells with all points on the surface if the
+    //- if >0 : preserve cells with all points on the surface if the
     //  resulting volume after snapping (by approximation) is larger than
     //  minVolCollapseRatio times old volume (i.e. not collapsed to flat cell).
     //  If <0 : delete always.
-    //minVolCollapseRatio 0.5;
+    //minVolCollapseRatio 0.1;
 
     // Advanced
 

applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C

@@ -168,6 +168,12 @@ int main(int argc, char *argv[])
         "add exposed internal faces to specified patch instead of to "
         "'oldInternalFaces'"
     );
+    argList::addOption
+    (
+        "resultTime",
+        "time",
+        "specify a time for the resulting mesh"
+    );
     #include "setRootCase.H"
     #include "createTime.H"
     runTime.functionObjects().off();
@@ -178,10 +184,12 @@ int main(int argc, char *argv[])
     #include "createNamedMesh.H"
 
 
-    const word oldInstance = mesh.pointsInstance();
-
     const word setName = args[1];
-    const bool overwrite = args.optionFound("overwrite");
+
+    word resultInstance = mesh.pointsInstance();
+    const bool specifiedInstance =
+        args.optionFound("overwrite")
+     || args.optionReadIfPresent("resultTime", resultInstance);
 
 
     Info<< "Reading cell set from " << setName << endl << endl;
@@ -347,13 +355,14 @@ int main(int argc, char *argv[])
     // Write mesh and fields to new time
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    if (!overwrite)
+    if (!specifiedInstance)
     {
         runTime++;
     }
     else
     {
-        subsetter.subMesh().setInstance(oldInstance);
+        runTime.setTime(instant(resultInstance), 0);
+        subsetter.subMesh().setInstance(runTime.timeName());
     }
 
     Info<< "Writing subsetted mesh and fields to time " << runTime.timeName()

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C

@@ -165,7 +165,7 @@ void Foam::GAMGAgglomeration::compactLevels(const label nCreatedLevels)
                     << setw(8) << setprecision(4) << totFaceCellRatio/totNprocs
                     << setw(8) << setprecision(4) << maxFaceCellRatio
                     << "        "
-                    << setw(8) << totNInt/totNprocs
+                    << setw(8) << scalar(totNInt)/totNprocs
                     << setw(8) << maxNInt
                     << "        "
                     << setw(8) << setprecision(4) << totRatio/totNprocs

src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.C

@@ -97,7 +97,6 @@ Foam::tmp<Foam::scalarField> Foam::polyMeshTools::faceSkewness
 {
     const labelList& own = mesh.faceOwner();
     const labelList& nei = mesh.faceNeighbour();
-    const faceList& fcs = mesh.faces();
     const polyBoundaryMesh& pbm = mesh.boundaryMesh();
 
     tmp<scalarField> tskew(new scalarField(mesh.nFaces()));
@@ -154,31 +153,16 @@ Foam::tmp<Foam::scalarField> Foam::polyMeshTools::faceSkewness
             {
                 label faceI = pp.start() + i;
 
-                vector Cpf = fCtrs[faceI] - cellCtrs[own[faceI]];
+                skew[faceI] = primitiveMeshTools::boundaryFaceSkewness
+                (
+                    mesh,
+                    p,
+                    fCtrs,
+                    fAreas,
 
-                vector normal = fAreas[faceI];
+                    faceI,
-                normal /= mag(normal) + VSMALL;
+                    cellCtrs[own[faceI]]
-                vector d = normal*(normal & Cpf);
+                );
-
-
-                // Skewness vector
-                vector sv =
-                    Cpf
-                  - ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + VSMALL))*d;
-                vector svHat = sv/(mag(sv) + VSMALL);
-
-                // Normalisation distance calculated as the approximate distance
-                // from the face centre to the edge of the face in the direction
-                // of the skewness
-                scalar fd = 0.4*mag(d) + VSMALL;
-                const face& f = fcs[faceI];
-                forAll(f, pi)
-                {
-                    fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[faceI])));
-                }
-
-                // Normalised skewness
-                skew[faceI] = mag(sv)/fd;
             }
         }
     }

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -64,6 +64,44 @@ Foam::scalar Foam::primitiveMeshTools::faceSkewness
     return mag(sv)/fd;
 }
 
+Foam::scalar Foam::primitiveMeshTools::boundaryFaceSkewness
+(
+    const primitiveMesh& mesh,
+    const pointField& p,
+    const vectorField& fCtrs,
+    const vectorField& fAreas,
+
+    const label faceI,
+    const point& ownCc
+)
+{
+    vector Cpf = fCtrs[faceI] - ownCc;
+
+    vector normal = fAreas[faceI];
+    normal /= mag(normal) + VSMALL;
+    vector d = normal*(normal & Cpf);
+
+
+    // Skewness vector
+    vector sv =
+        Cpf
+      - ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + VSMALL))*d;
+    vector svHat = sv/(mag(sv) + VSMALL);
+
+    // Normalisation distance calculated as the approximate distance
+    // from the face centre to the edge of the face in the direction
+    // of the skewness
+    scalar fd = 0.4*mag(d) + VSMALL;
+    const face& f = mesh.faces()[faceI];
+    forAll(f, pi)
+    {
+        fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[faceI])));
+    }
+
+    // Normalised skewness
+    return mag(sv)/fd;
+}
+
 
 Foam::scalar Foam::primitiveMeshTools::faceOrthogonality
 (
@@ -119,7 +157,6 @@ Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::faceSkewness
 {
     const labelList& own = mesh.faceOwner();
     const labelList& nei = mesh.faceNeighbour();
-    const faceList& fcs = mesh.faces();
 
     tmp<scalarField> tskew(new scalarField(mesh.nFaces()));
     scalarField& skew = tskew();
@@ -145,31 +182,15 @@ Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::faceSkewness
 
     for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
     {
-        vector Cpf = fCtrs[faceI] - cellCtrs[own[faceI]];
+        skew[faceI] = boundaryFaceSkewness
-
+        (
-        vector normal = fAreas[faceI];
+            mesh,
-        normal /= mag(normal) + VSMALL;
+            p,
-        vector d = normal*(normal & Cpf);
+            fCtrs,
-
+            fAreas,
-
+            faceI,
-        // Skewness vector
+            cellCtrs[own[faceI]]
-        vector sv =
+        );
-            Cpf
-          - ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + VSMALL))*d;
-        vector svHat = sv/(mag(sv) + VSMALL);
-
-        // Normalisation distance calculated as the approximate distance
-        // from the face centre to the edge of the face in the direction
-        // of the skewness
-        scalar fd = 0.4*mag(d) + VSMALL;
-        const face& f = fcs[faceI];
-        forAll(f, pi)
-        {
-            fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[faceI])));
-        }
-
-        // Normalised skewness
-        skew[faceI] = mag(sv)/fd;
     }
 
     return tskew;

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.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
@@ -48,32 +48,6 @@ namespace Foam
 
 class primitiveMeshTools
 {
-
-protected:
-
-   // Protected Member Functions
-
-        //- Skewness of single face
-        static scalar faceSkewness
-        (
-            const primitiveMesh& mesh,
-            const pointField& p,
-            const vectorField& fCtrs,
-            const vectorField& fAreas,
-
-            const label faceI,
-            const point& ownCc,
-            const point& neiCc
-        );
-
-        //- Orthogonality of single face
-        static scalar faceOrthogonality
-        (
-            const point& ownCc,
-            const point& neiCc,
-            const vector& s
-        );
-
 public:
 
     //- Generate non-orthogonality field (internal faces only)
@@ -154,6 +128,41 @@ public:
         const PackedBoolList& internalOrCoupledFace
     );
 
+
+    // Helpers: single face check
+
+        //- Skewness of single face
+        static scalar faceSkewness
+        (
+            const primitiveMesh& mesh,
+            const pointField& p,
+            const vectorField& fCtrs,
+            const vectorField& fAreas,
+
+            const label faceI,
+            const point& ownCc,
+            const point& neiCc
+        );
+
+        //- Skewness of single boundary face
+        static scalar boundaryFaceSkewness
+        (
+            const primitiveMesh& mesh,
+            const pointField& p,
+            const vectorField& fCtrs,
+            const vectorField& fAreas,
+
+            const label faceI,
+            const point& ownCc
+        );
+
+        //- Orthogonality of single face
+        static scalar faceOrthogonality
+        (
+            const point& ownCc,
+            const point& neiCc,
+            const vector& s
+        );
 };
 
 

src/dynamicMesh/motionSmoother/motionSmootherCheck.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
@@ -239,6 +239,7 @@ bool Foam::motionSmoother::checkMesh
             maxIntSkew,
             maxBounSkew,
             mesh,
+            mesh.points(),
             mesh.cellCentres(),
             mesh.faceCentres(),
             mesh.faceAreas(),
@@ -481,14 +482,14 @@ bool Foam::motionSmoother::checkMesh
     (
         readScalar(dict.lookup("minArea", true))
     );
-    const scalar maxIntSkew
+    //const scalar maxIntSkew
-    (
+    //(
-        readScalar(dict.lookup("maxInternalSkewness", true))
+    //    readScalar(dict.lookup("maxInternalSkewness", true))
-    );
+    //);
-    const scalar maxBounSkew
+    //const scalar maxBounSkew
-    (
+    //(
-        readScalar(dict.lookup("maxBoundarySkewness", true))
+    //    readScalar(dict.lookup("maxBoundarySkewness", true))
-    );
+    //);
     const scalar minWeight
     (
         readScalar(dict.lookup("minFaceWeight", true))
@@ -615,27 +616,30 @@ bool Foam::motionSmoother::checkMesh
         nWrongFaces = nNewWrongFaces;
     }
 
-    if (maxIntSkew > 0 || maxBounSkew > 0)
-    {
-        meshGeom.checkFaceSkewness
-        (
-            report,
-            maxIntSkew,
-            maxBounSkew,
-            checkFaces,
-            baffles,
-            &wrongFaces
-        );
 
-        label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
+    //- Note: cannot check the skewness without the points and don't want
-
+    //  to store them on polyMeshGeometry.
-        Info<< "    faces with skewness > "
+    //if (maxIntSkew > 0 || maxBounSkew > 0)
-            << setw(3) << maxIntSkew
+    //{
-            << " (internal) or " << setw(3) << maxBounSkew
+    //    meshGeom.checkFaceSkewness
-            << " (boundary) : " << nNewWrongFaces-nWrongFaces << endl;
+    //    (
-
+    //        report,
-        nWrongFaces = nNewWrongFaces;
+    //        maxIntSkew,
-    }
+    //        maxBounSkew,
+    //        checkFaces,
+    //        baffles,
+    //        &wrongFaces
+    //    );
+    //
+    //    label nNewWrongFaces = returnReduce(wrongFaces.size(),sumOp<label>());
+    //
+    //    Info<< "    faces with skewness > "
+    //        << setw(3) << maxIntSkew
+    //        << " (internal) or " << setw(3) << maxBounSkew
+    //        << " (boundary) : " << nNewWrongFaces-nWrongFaces << endl;
+    //
+    //    nWrongFaces = nNewWrongFaces;
+    //}
 
     if (minWeight >= 0 && minWeight < 1)
     {

src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C

@@ -29,6 +29,7 @@ License
 #include "tetrahedron.H"
 #include "syncTools.H"
 #include "unitConversion.H"
+#include "primitiveMeshTools.H"
 
 namespace Foam
 {
@@ -286,29 +287,6 @@ Foam::scalar Foam::polyMeshGeometry::checkNonOrtho
 }
 
 
-Foam::scalar Foam::polyMeshGeometry::calcSkewness
-(
-    const point& ownCc,
-    const point& neiCc,
-    const point& fc
-)
-{
-    scalar dOwn = mag(fc - ownCc);
-    scalar dNei = mag(fc - neiCc);
-
-    point faceIntersection =
-        ownCc*dNei/(dOwn+dNei+VSMALL)
-      + neiCc*dOwn/(dOwn+dNei+VSMALL);
-
-    return
-        mag(fc - faceIntersection)
-      / (
-            mag(neiCc-ownCc)
-          + VSMALL
-        );
-}
-
-
 // Create the neighbour pyramid - it will have positive volume
 bool Foam::polyMeshGeometry::checkFaceTet
 (
@@ -1008,6 +986,7 @@ bool Foam::polyMeshGeometry::checkFaceSkewness
     const scalar internalSkew,
     const scalar boundarySkew,
     const polyMesh& mesh,
+    const pointField& points,
     const vectorField& cellCentres,
     const vectorField& faceCentres,
     const vectorField& faceAreas,
@@ -1024,14 +1003,8 @@ bool Foam::polyMeshGeometry::checkFaceSkewness
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
     // Calculate coupled cell centre
-    pointField neiCc(mesh.nFaces()-mesh.nInternalFaces());
+    pointField neiCc;
-
+    syncTools::swapBoundaryCellPositions(mesh, cellCentres, neiCc);
-    for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
-    {
-        neiCc[faceI-mesh.nInternalFaces()] = cellCentres[own[faceI]];
-    }
-    syncTools::swapBoundaryFacePositions(mesh, neiCc);
-
 
     scalar maxSkew = 0;
 
@@ -1043,11 +1016,16 @@ bool Foam::polyMeshGeometry::checkFaceSkewness
 
         if (mesh.isInternalFace(faceI))
         {
-            scalar skewness = calcSkewness
+            scalar skewness = primitiveMeshTools::faceSkewness
             (
+                mesh,
+                points,
+                faceCentres,
+                faceAreas,
+
+                faceI,
                 cellCentres[own[faceI]],
-                cellCentres[nei[faceI]],
+                cellCentres[nei[faceI]]
-                faceCentres[faceI]
             );
 
             // Check if the skewness vector is greater than the PN vector.
@@ -1073,11 +1051,16 @@ bool Foam::polyMeshGeometry::checkFaceSkewness
         }
         else if (patches[patches.whichPatch(faceI)].coupled())
         {
-            scalar skewness = calcSkewness
+            scalar skewness = primitiveMeshTools::faceSkewness
             (
+                mesh,
+                points,
+                faceCentres,
+                faceAreas,
+
+                faceI,
                 cellCentres[own[faceI]],
-                neiCc[faceI-mesh.nInternalFaces()],
+                neiCc[faceI-mesh.nInternalFaces()]
-                faceCentres[faceI]
             );
 
             // Check if the skewness vector is greater than the PN vector.
@@ -1103,21 +1086,17 @@ bool Foam::polyMeshGeometry::checkFaceSkewness
         }
         else
         {
-            // Boundary faces: consider them to have only skewness error.
+            scalar skewness = primitiveMeshTools::boundaryFaceSkewness
-            // (i.e. treat as if mirror cell on other side)
+            (
+                mesh,
+                points,
+                faceCentres,
+                faceAreas,
 
-            vector faceNormal = faceAreas[faceI];
+                faceI,
-            faceNormal /= mag(faceNormal) + ROOTVSMALL;
+                cellCentres[own[faceI]]
+            );
 
-            vector dOwn = faceCentres[faceI] - cellCentres[own[faceI]];
-
-            vector dWall = faceNormal*(faceNormal & dOwn);
-
-            point faceIntersection = cellCentres[own[faceI]] + dWall;
-
-            scalar skewness =
-                mag(faceCentres[faceI] - faceIntersection)
-                /(2*mag(dWall) + ROOTVSMALL);
 
             // Check if the skewness vector is greater than the PN vector.
             // This does not cause trouble but is a good indication of a poor
@@ -1148,12 +1127,18 @@ bool Foam::polyMeshGeometry::checkFaceSkewness
         label face1 = baffles[i].second();
 
         const point& ownCc = cellCentres[own[face0]];
+        const point& neiCc = cellCentres[own[face1]];
 
-        scalar skewness = calcSkewness
+        scalar skewness = primitiveMeshTools::faceSkewness
         (
+            mesh,
+            points,
+            faceCentres,
+            faceAreas,
+
+            face0,
             ownCc,
-            cellCentres[own[face1]],
+            neiCc
-            faceCentres[face0]
         );
 
         // Check if the skewness vector is greater than the PN vector.
@@ -2361,30 +2346,30 @@ bool Foam::polyMeshGeometry::checkFaceTets
 }
 
 
-bool Foam::polyMeshGeometry::checkFaceSkewness
+//bool Foam::polyMeshGeometry::checkFaceSkewness
-(
+//(
-    const bool report,
+//    const bool report,
-    const scalar internalSkew,
+//    const scalar internalSkew,
-    const scalar boundarySkew,
+//    const scalar boundarySkew,
-    const labelList& checkFaces,
+//    const labelList& checkFaces,
-    const List<labelPair>& baffles,
+//    const List<labelPair>& baffles,
-    labelHashSet* setPtr
+//    labelHashSet* setPtr
-) const
+//) const
-{
+//{
-    return checkFaceSkewness
+//    return checkFaceSkewness
-    (
+//    (
-        report,
+//        report,
-        internalSkew,
+//        internalSkew,
-        boundarySkew,
+//        boundarySkew,
-        mesh_,
+//        mesh_,
-        cellCentres_,
+//        cellCentres_,
-        faceCentres_,
+//        faceCentres_,
-        faceAreas_,
+//        faceAreas_,
-        checkFaces,
+//        checkFaces,
-        baffles,
+//        baffles,
-        setPtr
+//        setPtr
-    );
+//    );
-}
+//}
 
 
 bool Foam::polyMeshGeometry::checkFaceWeights

src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.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
@@ -229,6 +229,7 @@ public:
                 const scalar internalSkew,
                 const scalar boundarySkew,
                 const polyMesh& mesh,
+                const pointField& points,
                 const vectorField& cellCentres,
                 const vectorField& faceCentres,
                 const vectorField& faceAreas,
@@ -372,16 +373,6 @@ public:
                 labelHashSet* setPtr
             ) const;
 
-            bool checkFaceSkewness
-            (
-                const bool report,
-                const scalar internalSkew,
-                const scalar boundarySkew,
-                const labelList& checkFaces,
-                const List<labelPair>& baffles,
-                labelHashSet* setPtr
-            ) const;
-
             bool checkFaceWeights
             (
                 const bool report,

src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.C

@@ -130,8 +130,8 @@ void Foam::porosityModels::DarcyForchheimer::calcForce
     vectorField& force
 ) const
 {
-    scalarField Udiag(U.size());
+    scalarField Udiag(U.size(), 0.0);
-    vectorField Usource(U.size());
+    vectorField Usource(U.size(), vector::zero);
     const scalarField& V = mesh_.V();
 
     apply(Udiag, Usource, V, rho, mu, U);

src/finiteVolume/cfdTools/general/porosityModel/fixedCoeff/fixedCoeff.C

@@ -183,8 +183,8 @@ void Foam::porosityModels::fixedCoeff::calcForce
     vectorField& force
 ) const
 {
-    scalarField Udiag(U.size());
+    scalarField Udiag(U.size(), 0.0);
-    vectorField Usource(U.size());
+    vectorField Usource(U.size(), vector::zero);
     const scalarField& V = mesh_.V();
     scalar rhoRef = readScalar(coeffs_.lookup("rhoRef"));
 

src/finiteVolume/cfdTools/general/porosityModel/powerLaw/powerLaw.C

@@ -74,7 +74,7 @@ void Foam::porosityModels::powerLaw::calcForce
     vectorField& force
 ) const
 {
-    scalarField Udiag(U.size());
+    scalarField Udiag(U.size(), 0.0);
     const scalarField& V = mesh_.V();
 
     apply(Udiag, V, rho, U);

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
 {

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
+            //- Manipulate matrix
-            virtual void updateInterfaceMatrix
+            virtual void manipulateMatrix(fvMatrix<Type>& matrix);
-            (
-                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;
 
 
         // Cyclic AMI coupled interface functions

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;
 }
 
 

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()
+            virtual void updateCoeffs();
-            {
+
-                updated_ = true;
+            //- 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
 

src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceMDLimitedGrad/faceMDLimitedGrads.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -235,7 +235,7 @@ Foam::fv::faceMDLimitedGrad<Foam::vector>::calcGrad
             g[own],
             maxFace - vvfOwn,
             minFace - vvfOwn,
-            Cf[facei] - C[nei]
+            Cf[facei] - C[own]
         );
 
 

src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwind.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -111,16 +111,7 @@ public:
                     mesh.gradScheme(gradSchemeName_)
                 )
             )
-        {
+        {}
-            if (!schemeData.eof())
-            {
-                IOWarningIn("linearUpwind(const fvMesh&, Istream&)", schemeData)
-                    << "unexpected additional entries in stream." << nl
-                    << "    Only the name of the gradient scheme in the"
-                       " 'gradSchemes' dictionary should be specified."
-                    << endl;
-            }
-        }
 
         //- Construct from faceFlux and Istream
         linearUpwind
@@ -140,21 +131,7 @@ public:
                     mesh.gradScheme(gradSchemeName_)
                 )
             )
-        {
+        {}
-
-            if (!schemeData.eof())
-            {
-                IOWarningIn
-                (
-                    "linearUpwind(const fvMesh&, "
-                    "const surfaceScalarField& faceFlux, Istream&)",
-                    schemeData
-                )   << "unexpected additional entries in stream." << nl
-                    << "    Only the name of the gradient scheme in the"
-                       " 'gradSchemes' dictionary should be specified."
-                    << endl;
-            }
-        }
 
 
     // Member Functions

src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwindV.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -110,19 +110,7 @@ public:
                     mesh.gradScheme(gradSchemeName_)
                 )
             )
-        {
+        {}
-            if (!schemeData.eof())
-            {
-                IOWarningIn
-                (
-                    "linearUpwindV(const fvMesh&, Istream&)",
-                    schemeData
-                )   << "unexpected additional entries in stream." << nl
-                    << "    Only the name of the gradient scheme in the"
-                       " 'gradSchemes' dictionary should be specified."
-                    << endl;
-            }
-        }
 
         //- Construct from faceFlux and Istream
         linearUpwindV
@@ -142,20 +130,7 @@ public:
                     mesh.gradScheme(gradSchemeName_)
                 )
             )
-        {
+        {}
-            if (!schemeData.eof())
-            {
-                IOWarningIn
-                (
-                    "linearUpwindV(const fvMesh&, "
-                    "const surfaceScalarField& faceFlux, Istream&)",
-                    schemeData
-                )   << "unexpected additional entries in stream." << nl
-                    << "    Only the name of the gradient scheme in the"
-                       " 'gradSchemes' dictionary should be specified."
-                    << endl;
-            }
-        }
 
 
     // Member Functions

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

@@ -336,7 +336,7 @@ void Foam::autoLayerDriver::smoothNormals
 ) const
 {
     // Get smoothly varying internal normals field.
-    Info<< "shrinkMeshDistance : Smoothing normals ..." << endl;
+    Info<< "shrinkMeshDistance : Smoothing normals in interior ..." << endl;
 
     const fvMesh& mesh = meshRefiner_.mesh();
     const edgeList& edges = mesh.edges();
@@ -1161,6 +1161,11 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
         scalar mDist = medialDist[pointI];
 
         if (wDist2 < sqr(SMALL) && mDist < SMALL)
+        //- Note: maybe less strict:
+        //(
+        //    wDist2 < sqr(meshRefiner_.mergeDistance())
+        // && mDist < meshRefiner_.mergeDistance()
+        //)
         {
             medialRatio[pointI] = 0.0;
         }

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

@@ -35,6 +35,7 @@ License
 #include "shellSurfaces.H"
 #include "mapDistributePolyMesh.H"
 #include "unitConversion.H"
+#include "snapParameters.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -610,9 +611,16 @@ Foam::label Foam::autoRefineDriver::danglingCellRefine
             << " cells (out of " << mesh.globalData().nTotalCells()
             << ')' << endl;
 
-        // Stop when no cells to refine. No checking of minRefineCells since
+        // Stop when no cells to refine. After a few iterations check if too
-        // too few cells
+        // few cells
-        if (nCellsToRefine == 0)
+        if
+        (
+            nCellsToRefine == 0
+         || (
+                iter >= 1
+             && nCellsToRefine <= refineParams.minRefineCells()
+            )
+        )
         {
             Info<< "Stopping refining since too few cells selected."
                 << nl << endl;
@@ -870,6 +878,7 @@ Foam::label Foam::autoRefineDriver::shellRefine
 void Foam::autoRefineDriver::baffleAndSplitMesh
 (
     const refinementParameters& refineParams,
+    const snapParameters& snapParams,
     const bool handleSnapProblems,
     const dictionary& motionDict
 )
@@ -887,10 +896,17 @@ void Foam::autoRefineDriver::baffleAndSplitMesh
     meshRefiner_.baffleAndSplitMesh
     (
         handleSnapProblems,             // detect&remove potential snap problem
+
+        // Snap problem cell detection
+        snapParams,
+        refineParams.useTopologicalSnapDetection(),
         false,                          // perpendicular edge connected cells
         scalarField(0),                 // per region perpendicular angle
+
+        // Free standing baffles
         !handleSnapProblems,            // merge free standing baffles?
         refineParams.planarAngle(),
+
         motionDict,
         const_cast<Time&>(mesh.time()),
         globalToMasterPatch_,
@@ -951,6 +967,7 @@ void Foam::autoRefineDriver::zonify
 void Foam::autoRefineDriver::splitAndMergeBaffles
 (
     const refinementParameters& refineParams,
+    const snapParameters& snapParams,
     const bool handleSnapProblems,
     const dictionary& motionDict
 )
@@ -973,10 +990,17 @@ void Foam::autoRefineDriver::splitAndMergeBaffles
     meshRefiner_.baffleAndSplitMesh
     (
         handleSnapProblems,
+
+        // Snap problem cell detection
+        snapParams,
+        refineParams.useTopologicalSnapDetection(),
         handleSnapProblems,                 // remove perp edge connected cells
         perpAngle,                          // perp angle
+
+        // Free standing baffles
         true,                               // merge free standing baffles?
         refineParams.planarAngle(),         // planar angle
+
         motionDict,
         const_cast<Time&>(mesh.time()),
         globalToMasterPatch_,
@@ -1081,6 +1105,7 @@ void Foam::autoRefineDriver::doRefine
 (
     const dictionary& refineDict,
     const refinementParameters& refineParams,
+    const snapParameters& snapParams,
     const bool prepareForSnapping,
     const dictionary& motionDict
 )
@@ -1146,13 +1171,25 @@ void Foam::autoRefineDriver::doRefine
 
     // Introduce baffles at surface intersections. Remove sections unreachable
     // from keepPoint.
-    baffleAndSplitMesh(refineParams, prepareForSnapping, motionDict);
+    baffleAndSplitMesh
+    (
+        refineParams,
+        snapParams,
+        prepareForSnapping,
+        motionDict
+    );
 
     // Mesh is at its finest. Do optional zoning.
     zonify(refineParams);
 
     // Pull baffles apart
-    splitAndMergeBaffles(refineParams, prepareForSnapping, motionDict);
+    splitAndMergeBaffles
+    (
+        refineParams,
+        snapParams,
+        prepareForSnapping,
+        motionDict
+    );
 
     // Do something about cells with refined faces on the boundary
     if (prepareForSnapping)

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

@@ -43,6 +43,8 @@ namespace Foam
 
 // Forward declaration of classes
 class refinementParameters;
+class snapParameters;
+
 class meshRefinement;
 class decompositionMethod;
 class fvMeshDistribute;
@@ -121,6 +123,7 @@ class autoRefineDriver
         void baffleAndSplitMesh
         (
             const refinementParameters& refineParams,
+            const snapParameters& snapParams,
             const bool handleSnapProblems,
             const dictionary& motionDict
         );
@@ -131,6 +134,7 @@ class autoRefineDriver
         void splitAndMergeBaffles
         (
             const refinementParameters& refineParams,
+            const snapParameters& snapParams,
             const bool handleSnapProblems,
             const dictionary& motionDict
         );
@@ -176,6 +180,7 @@ public:
         (
             const dictionary& refineDict,
             const refinementParameters& refineParams,
+            const snapParameters& snapParams,
             const bool prepareForSnapping,
             const dictionary& motionDict
         );

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

@@ -121,7 +121,7 @@ Foam::pointField Foam::autoSnapDriver::smoothPatchDisplacement
 (
     const motionSmoother& meshMover,
     const List<labelPair>& baffles
-) const
+)
 {
     const indirectPrimitivePatch& pp = meshMover.patch();
 
@@ -615,14 +615,14 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::mergeZoneBaffles
 
 Foam::scalarField Foam::autoSnapDriver::calcSnapDistance
 (
+    const fvMesh& mesh,
     const snapParameters& snapParams,
     const indirectPrimitivePatch& pp
-) const
+)
 {
     const edgeList& edges = pp.edges();
     const labelListList& pointEdges = pp.pointEdges();
     const pointField& localPoints = pp.localPoints();
-    const fvMesh& mesh = meshRefiner_.mesh();
 
     scalarField maxEdgeLen(localPoints.size(), -GREAT);
 
@@ -655,13 +655,14 @@ Foam::scalarField Foam::autoSnapDriver::calcSnapDistance
 
 void Foam::autoSnapDriver::preSmoothPatch
 (
+    const meshRefinement& meshRefiner,
     const snapParameters& snapParams,
     const label nInitErrors,
     const List<labelPair>& baffles,
     motionSmoother& meshMover
-) const
+)
 {
-    const fvMesh& mesh = meshRefiner_.mesh();
+    const fvMesh& mesh = meshRefiner.mesh();
 
     labelList checkFaces;
 
@@ -724,11 +725,11 @@ void Foam::autoSnapDriver::preSmoothPatch
     {
         const_cast<Time&>(mesh.time())++;
         Info<< "Writing patch smoothed mesh to time "
-            << meshRefiner_.timeName() << '.' << endl;
+            << meshRefiner.timeName() << '.' << endl;
-        meshRefiner_.write
+        meshRefiner.write
         (
             debug,
-            mesh.time().path()/meshRefiner_.timeName()
+            mesh.time().path()/meshRefiner.timeName()
         );
         Info<< "Dumped mesh in = "
             << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
@@ -742,12 +743,11 @@ void Foam::autoSnapDriver::preSmoothPatch
 // Get (pp-local) indices of points that are both on zone and on patched surface
 Foam::labelList Foam::autoSnapDriver::getZoneSurfacePoints
 (
+    const fvMesh& mesh,
     const indirectPrimitivePatch& pp,
     const word& zoneName
-) const
+)
 {
-    const fvMesh& mesh = meshRefiner_.mesh();
-
     label zoneI = mesh.faceZones().findZoneID(zoneName);
 
     if (zoneI == -1)
@@ -755,7 +755,7 @@ Foam::labelList Foam::autoSnapDriver::getZoneSurfacePoints
         FatalErrorIn
         (
             "autoSnapDriver::getZoneSurfacePoints"
-            "(const indirectPrimitivePatch&, const word&)"
+            "(const fvMesh&, const indirectPrimitivePatch&, const word&)"
         )   << "Cannot find zone " << zoneName
             << exit(FatalError);
     }
@@ -793,17 +793,17 @@ Foam::labelList Foam::autoSnapDriver::getZoneSurfacePoints
 
 Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
 (
+    const meshRefinement& meshRefiner,
     const scalarField& snapDist,
-    motionSmoother& meshMover
+    const indirectPrimitivePatch& pp
-) const
+)
 {
     Info<< "Calculating patchDisplacement as distance to nearest surface"
         << " point ..." << endl;
 
-    const indirectPrimitivePatch& pp = meshMover.patch();
     const pointField& localPoints = pp.localPoints();
-    const refinementSurfaces& surfaces = meshRefiner_.surfaces();
+    const refinementSurfaces& surfaces = meshRefiner.surfaces();
-    const fvMesh& mesh = meshRefiner_.mesh();
+    const fvMesh& mesh = meshRefiner.mesh();
 
     // Displacement per patch point
     vectorField patchDisp(localPoints.size(), vector::zero);
@@ -815,9 +815,9 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
 
         // Divide surfaces into zoned and unzoned
         labelList zonedSurfaces =
-            meshRefiner_.surfaces().getNamedSurfaces();
+            meshRefiner.surfaces().getNamedSurfaces();
         labelList unzonedSurfaces =
-            meshRefiner_.surfaces().getUnnamedSurfaces();
+            meshRefiner.surfaces().getUnnamedSurfaces();
 
 
         // 1. All points to non-interface surfaces
@@ -870,6 +870,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
             (
                 getZoneSurfacePoints
                 (
+                    mesh,
                     pp,
                     faceZoneNames[zoneSurfI]
                 )
@@ -966,6 +967,254 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
 }
 
 
+////XXXXXXXXX
+//// Get (pp-local) indices of points that are on both patches
+//Foam::labelList Foam::autoSnapDriver::getPatchSurfacePoints
+//(
+//    const fvMesh& mesh,
+//    const indirectPrimitivePatch& allPp,
+//    const polyPatch& pp
+//)
+//{
+//    // Could use PrimitivePatch & localFaces to extract points but might just
+//    // as well do it ourselves.
+//
+//    boolList pointOnZone(allPp.nPoints(), false);
+//
+//    forAll(pp, i)
+//    {
+//        const face& f = pp[i];
+//
+//        forAll(f, fp)
+//        {
+//            label meshPointI = f[fp];
+//
+//            Map<label>::const_iterator iter =
+//                allPp.meshPointMap().find(meshPointI);
+//
+//            if (iter != allPp.meshPointMap().end())
+//            {
+//                label pointI = iter();
+//                pointOnZone[pointI] = true;
+//            }
+//        }
+//    }
+//
+//    return findIndices(pointOnZone, true);
+//}
+//Foam::vectorField Foam::autoSnapDriver::calcNearestLocalSurface
+//(
+//    const meshRefinement& meshRefiner,
+//    const scalarField& snapDist,
+//    const indirectPrimitivePatch& pp
+//)
+//{
+//    Info<< "Calculating patchDisplacement as distance to nearest"
+//        << " local surface point ..." << endl;
+//
+//    const pointField& localPoints = pp.localPoints();
+//    const refinementSurfaces& surfaces = meshRefiner.surfaces();
+//    const fvMesh& mesh = meshRefiner.mesh();
+//    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+//
+//
+////    // Assume that all patch-internal points get attracted to their surface
+////    // only. So we want to know if point is on multiple regions
+////
+////    labelList minPatch(mesh.nPoints(), labelMax);
+////    labelList maxPatch(mesh.nPoints(), labelMin);
+////
+////    forAll(meshMover.adaptPatchIDs(), i)
+////    {
+////        label patchI = meshMover.adaptPatchIDs()[i];
+////        const labelList& meshPoints = pbm[patchI].meshPoints();
+////
+////        forAll(meshPoints, meshPointI)
+////        {
+////            label meshPointI = meshPoints[meshPointI];
+////            minPatch[meshPointI] = min(minPatch[meshPointI], patchI);
+////            maxPatch[meshPointI] = max(maxPatch[meshPointI], patchI);
+////        }
+////    }
+////
+////    syncTools::syncPointList
+////    (
+////        mesh,
+////        minPatch,
+////        minEqOp<label>(),   // combine op
+////        labelMax            // null value
+////    );
+////    syncTools::syncPointList
+////    (
+////        mesh,
+////        maxPatch,
+////        maxEqOp<label>(),   // combine op
+////        labelMin            // null value
+////    );
+//
+//    // Now all points with minPatch != maxPatch will be on the outside of
+//    // the patch.
+//
+//    // Displacement per patch point
+//    vectorField patchDisp(localPoints.size(), vector::zero);
+//    // Current best snap distance
+//    scalarField minSnapDist(snapDist);
+//    // Current surface snapped to
+//    labelList snapSurf(localPoints.size(), -1);
+//
+//    const labelList& surfaceGeometry = surfaces.surfaces();
+//    forAll(surfaceGeometry, surfI)
+//    {
+//        label geomI = surfaceGeometry[surfI];
+//        const wordList& regNames = allGeometry.regionNames()[geomI];
+//        forAll(regNames, regionI)
+//        {
+//            label globalRegionI = surfaces.globalRegion(surfI, regionI);
+//            // Collect master patch points
+//            label masterPatchI = globalToMasterPatch_[globalRegionI];
+//            label slavePatchI = globalToSlavePatch_[globalRegionI];
+//
+//            labelList patchPointIndices
+//            (
+//                getPatchSurfacePoints
+//                (
+//                    mesh,
+//                    pp,
+//                    pbm[masterPatchI]
+//                )
+//            );
+//
+//            // Find nearest for points both on faceZone and pp.
+//            List<pointIndexHit> hitInfo;
+//            surfaces.findNearest
+//            (
+//                surfI,
+//                regionI,
+//                pointField(localPoints, patchPointIndices),
+//                sqr(scalarField(minSnapDist, patchPointIndices)),
+//                hitInfo
+//            );
+//
+//            forAll(hitInfo, i)
+//            {
+//                label pointI = patchPointIndices[i];
+//
+//                if (hitInfo[i].hit())
+//                {
+//                    const point& pt = hitInfo[i].hitPoint();
+//                    patchDisp[pointI] = pt-localPoints[pointI];
+//                    minSnapDist[pointI] = min
+//                    (
+//                        minSnapDist[pointI],
+//                        mag(patchDisp[pointI])
+//                    );
+//                    snapSurf[pointI] = surfI;
+//                }
+//            }
+//
+//            // Slave patch
+//            if (slavePatchI != masterPatchI)
+//            {
+//                labelList patchPointIndices
+//                (
+//                    getPatchSurfacePoints
+//                    (
+//                        mesh,
+//                        pp,
+//                        pbm[slavePatchI]
+//                    )
+//                );
+//
+//                // Find nearest for points both on faceZone and pp.
+//                List<pointIndexHit> hitInfo;
+//                surfaces.findNearest
+//                (
+//                    surfI,
+//                    regionI,
+//                    pointField(localPoints, patchPointIndices),
+//                    sqr(scalarField(minSnapDist, patchPointIndices)),
+//                    hitInfo
+//                );
+//
+//                forAll(hitInfo, i)
+//                {
+//                    label pointI = patchPointIndices[i];
+//
+//                    if (hitInfo[i].hit())
+//                    {
+//                        const point& pt = hitInfo[i].hitPoint();
+//                        patchDisp[pointI] = pt-localPoints[pointI];
+//                        minSnapDist[pointI] = min
+//                        (
+//                            minSnapDist[pointI],
+//                            mag(patchDisp[pointI])
+//                        );
+//                        snapSurf[pointI] = surfI;
+//                    }
+//                }
+//            }
+//        }
+//    }
+//
+//
+//    // Check if all points are being snapped
+//    forAll(snapSurf, pointI)
+//    {
+//        if (snapSurf[pointI] == -1)
+//        {
+//            WarningIn("autoSnapDriver::calcNearestLocalSurface(..)")
+//                << "For point:" << pointI
+//                << " coordinate:" << localPoints[pointI]
+//                << " did not find any surface within:"
+//                << minSnapDist[pointI]
+//                << " metre." << endl;
+//        }
+//    }
+//
+//    {
+//        scalarField magDisp(mag(patchDisp));
+//
+//        Info<< "Wanted displacement : average:"
+//            << gSum(magDisp)/returnReduce(patchDisp.size(), sumOp<label>())
+//            << " min:" << gMin(magDisp)
+//            << " max:" << gMax(magDisp) << endl;
+//    }
+//
+//
+//    Info<< "Calculated surface displacement in = "
+//        << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
+//
+//
+//    // Limit amount of movement.
+//    forAll(patchDisp, patchPointI)
+//    {
+//        scalar magDisp = mag(patchDisp[patchPointI]);
+//
+//        if (magDisp > snapDist[patchPointI])
+//        {
+//            patchDisp[patchPointI] *= snapDist[patchPointI] / magDisp;
+//
+//            Pout<< "Limiting displacement for " << patchPointI
+//                << " from " << magDisp << " to " << snapDist[patchPointI]
+//                << endl;
+//        }
+//    }
+//
+//    // Points on zones in one domain but only present as point on other
+//    // will not do condition 2 on all. Sync explicitly.
+//    syncTools::syncPointList
+//    (
+//        mesh,
+//        pp.meshPoints(),
+//        patchDisp,
+//        minMagSqrEqOp<point>(),    // combine op
+//        vector(GREAT, GREAT, GREAT)// null value (note: cannot use VGREAT)
+//    );
+//
+//    return patchDisp;
+//}
+////XXXXXXXXX
+
 void Foam::autoSnapDriver::smoothDisplacement
 (
     const snapParameters& snapParams,
@@ -1153,7 +1402,15 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::repatchToSurface
         scalarField faceSnapDist(pp.size(), -GREAT);
         {
             // Distance to attract to nearest feature on surface
-            const scalarField snapDist(calcSnapDistance(snapParams, pp));
+            const scalarField snapDist
+            (
+                calcSnapDistance
+                (
+                    mesh,
+                    snapParams,
+                    pp
+                )
+            );
 
             const faceList& localFaces = pp.localFaces();
 
@@ -1429,7 +1686,7 @@ void Foam::autoSnapDriver::doSnap
         indirectPrimitivePatch& pp = ppPtr();
 
         // Distance to attract to nearest feature on surface
-        const scalarField snapDist(calcSnapDistance(snapParams, pp));
+        const scalarField snapDist(calcSnapDistance(mesh, snapParams, pp));
 
 
         // Construct iterative mesh mover.
@@ -1467,7 +1724,14 @@ void Foam::autoSnapDriver::doSnap
             << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
 
         // Pre-smooth patch vertices (so before determining nearest)
-        preSmoothPatch(snapParams, nInitErrors, baffles, meshMover);
+        preSmoothPatch
+        (
+            meshRefiner_,
+            snapParams,
+            nInitErrors,
+            baffles,
+            meshMover
+        );
 
 
         for (label iter = 0; iter < nFeatIter; iter++)
@@ -1478,7 +1742,12 @@ void Foam::autoSnapDriver::doSnap
 
             // Calculate displacement at every patch point. Insert into
             // meshMover.
-            vectorField disp = calcNearestSurface(snapDist, meshMover);
+            vectorField disp = calcNearestSurface
+            (
+                meshRefiner_,
+                snapDist,
+                meshMover.patch()
+            );
 
             // Override displacement with feature edge attempt
             if (doFeatures)

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -81,11 +81,11 @@ class autoSnapDriver
 
             //- Calculate displacement per patch point to smooth out patch.
             //  Quite complicated in determining which points to move where.
-            pointField smoothPatchDisplacement
+            static pointField smoothPatchDisplacement
             (
                 const motionSmoother&,
                 const List<labelPair>&
-            ) const;
+            );
 
             //- Check that face zones are synced
             void checkCoupledFaceZones() const;
@@ -406,37 +406,51 @@ public:
             autoPtr<mapPolyMesh> mergeZoneBaffles(const List<labelPair>&);
 
             //- Calculate edge length per patch point.
-            scalarField calcSnapDistance
+            static scalarField calcSnapDistance
             (
+                const fvMesh& mesh,
                 const snapParameters& snapParams,
                 const indirectPrimitivePatch&
-            ) const;
+            );
 
             //- Smooth the mesh (patch and internal) to increase visibility
             //  of surface points (on castellated mesh) w.r.t. surface.
-            void preSmoothPatch
+            static void preSmoothPatch
             (
+                const meshRefinement& meshRefiner,
                 const snapParameters& snapParams,
                 const label nInitErrors,
                 const List<labelPair>& baffles,
                 motionSmoother&
-            ) const;
+            );
 
             //- Get points both on patch and facezone.
-            labelList getZoneSurfacePoints
+            static labelList getZoneSurfacePoints
             (
+                const fvMesh& mesh,
                 const indirectPrimitivePatch&,
                 const word& zoneName
-            ) const;
+            );
 
             //- Per patch point calculate point on nearest surface. Set as
             //  boundary conditions of motionSmoother displacement field. Return
             //  displacement of patch points.
-            vectorField calcNearestSurface
+            static vectorField calcNearestSurface
             (
+                const meshRefinement& meshRefiner,
                 const scalarField& snapDist,
-                motionSmoother& meshMover
+                const indirectPrimitivePatch&
-            ) const;
+            );
+
+            ////- Per patch point calculate point on nearest surface. Set as
+            ////  boundary conditions of motionSmoother displacement field.
+            ////  Return displacement of patch points.
+            //static vectorField calcNearestLocalSurface
+            //(
+            //    const meshRefinement& meshRefiner,
+            //    const scalarField& snapDist,
+            //    const indirectPrimitivePatch&
+            //);
 
             //- Smooth the displacement field to the internal.
             void smoothDisplacement

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

@@ -45,7 +45,11 @@ Foam::refinementParameters::refinementParameters
     nBufferLayers_(readLabel(dict.lookup("nBufferLayers"))),
     keepPoints_(dict.lookup("keepPoints")),
     allowFreeStandingZoneFaces_(dict.lookup("allowFreeStandingZoneFaces")),
-    maxLoadUnbalance_(dict.lookupOrDefault<scalar>("maxLoadUnbalance",0))
+    useTopologicalSnapDetection_
+    (
+        dict.lookupOrDefault<bool>("useTopologicalSnapDetection", true)
+    ),
+    maxLoadUnbalance_(dict.lookupOrDefault<scalar>("maxLoadUnbalance", 0))
 {}
 
 
@@ -65,7 +69,11 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
     nBufferLayers_(readLabel(dict.lookup("nCellsBetweenLevels"))),
     keepPoints_(pointField(1, dict.lookup("locationInMesh"))),
     allowFreeStandingZoneFaces_(dict.lookup("allowFreeStandingZoneFaces")),
-    maxLoadUnbalance_(dict.lookupOrDefault<scalar>("maxLoadUnbalance",0))
+    useTopologicalSnapDetection_
+    (
+        dict.lookupOrDefault<bool>("useTopologicalSnapDetection", true)
+    ),
+    maxLoadUnbalance_(dict.lookupOrDefault<scalar>("maxLoadUnbalance", 0))
 {
     scalar featAngle(readScalar(dict.lookup("resolveFeatureAngle")));
 
@@ -83,7 +91,7 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::labelList Foam::refinementParameters::findCells(const polyMesh& mesh)
- const
+const
 {
     // Force calculation of tet-diag decomposition (for use in findCell)
     (void)mesh.tetBasePtIs();

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

@@ -80,6 +80,10 @@ class refinementParameters
         //  cellZone?
         Switch allowFreeStandingZoneFaces_;
 
+        //- Use old topology based problem-cell removal (cells with 8 points
+        //  on surface)
+        Switch useTopologicalSnapDetection_;
+
         //- Allowed load unbalance
         scalar maxLoadUnbalance_;
 
@@ -157,6 +161,13 @@ public:
                 return allowFreeStandingZoneFaces_;
             }
 
+            //- Use old topology based problem-cell removal
+            //  (cells with 8 points on surface)
+            bool useTopologicalSnapDetection() const
+            {
+                return useTopologicalSnapDetection_;
+            }
+
             //- Allowed load unbalance
             scalar maxLoadUnbalance() const
             {

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

@@ -72,6 +72,8 @@ class globalIndex;
 class removePoints;
 class localPointRegion;
 
+class snapParameters;
+
 /*---------------------------------------------------------------------------*\
                            Class meshRefinement Declaration
 \*---------------------------------------------------------------------------*/
@@ -470,6 +472,14 @@ private:
                 const labelList& globalToMasterPatch
             ) const;
 
+            //- Returns list with for every internal face -1 or the patch
+            //  they should be baffled into.
+            labelList markFacesOnProblemCellsGeometric
+            (
+                const snapParameters& snapParams,
+                const dictionary& motionDict
+            ) const;
+
 
         // Baffle merging
 
@@ -536,6 +546,8 @@ private:
             //- Remove any loose standing cells
             void handleSnapProblems
             (
+                const snapParameters& snapParams,
+                const bool useTopologicalSnapDetection,
                 const bool removeEdgeConnectedCells,
                 const scalarField& perpendicularAngle,
                 const dictionary& motionDict,
@@ -760,10 +772,17 @@ public:
             void baffleAndSplitMesh
             (
                 const bool handleSnapProblems,
+
+                // How to remove problem snaps
+                const snapParameters& snapParams,
+                const bool useTopologicalSnapDetection,
                 const bool removeEdgeConnectedCells,
                 const scalarField& perpendicularAngle,
+
+                // How to handle free-standing baffles
                 const bool mergeFreeStanding,
                 const scalar freeStandingAngle,
+
                 const dictionary& motionDict,
                 Time& runTime,
                 const labelList& globalToMasterPatch,

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

@@ -1872,6 +1872,8 @@ void Foam::meshRefinement::makeConsistentFaceIndex
 
 void Foam::meshRefinement::handleSnapProblems
 (
+    const snapParameters& snapParams,
+    const bool useTopologicalSnapDetection,
     const bool removeEdgeConnectedCells,
     const scalarField& perpendicularAngle,
     const dictionary& motionDict,
@@ -1885,44 +1887,39 @@ void Foam::meshRefinement::handleSnapProblems
         << "----------------------------------------------" << nl
         << endl;
 
-    labelList facePatch
+    labelList facePatch;
-    (
+    if (useTopologicalSnapDetection)
-        markFacesOnProblemCells
+    {
+        facePatch = markFacesOnProblemCells
         (
             motionDict,
             removeEdgeConnectedCells,
             perpendicularAngle,
             globalToMasterPatch
-        )
         );
+    }
+    else
+    {
+        facePatch = markFacesOnProblemCellsGeometric(snapParams, motionDict);
+    }
     Info<< "Analyzed problem cells in = "
         << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
 
     if (debug&meshRefinement::MESH)
     {
-        faceSet problemTopo(mesh_, "problemFacesTopo", 100);
+        faceSet problemFaces(mesh_, "problemFaces", 100);
 
-        const labelList facePatchTopo
+        forAll(facePatch, faceI)
-        (
-            markFacesOnProblemCells
-            (
-                motionDict,
-                removeEdgeConnectedCells,
-                perpendicularAngle,
-                globalToMasterPatch
-            )
-        );
-        forAll(facePatchTopo, faceI)
         {
-            if (facePatchTopo[faceI] != -1)
+            if (facePatch[faceI] != -1)
             {
-                problemTopo.insert(faceI);
+                problemFaces.insert(faceI);
             }
         }
-        problemTopo.instance() = timeName();
+        problemFaces.instance() = timeName();
-        Pout<< "Dumping " << problemTopo.size()
+        Pout<< "Dumping " << problemFaces.size()
-            << " problem faces to " << problemTopo.objectPath() << endl;
+            << " problem faces to " << problemFaces.objectPath() << endl;
-        problemTopo.write();
+        problemFaces.write();
     }
 
     Info<< "Introducing baffles to delete problem cells." << nl << endl;
@@ -1962,6 +1959,8 @@ void Foam::meshRefinement::handleSnapProblems
 void Foam::meshRefinement::baffleAndSplitMesh
 (
     const bool doHandleSnapProblems,
+    const snapParameters& snapParams,
+    const bool useTopologicalSnapDetection,
     const bool removeEdgeConnectedCells,
     const scalarField& perpendicularAngle,
     const bool mergeFreeStanding,
@@ -2029,83 +2028,10 @@ void Foam::meshRefinement::baffleAndSplitMesh
 
     if (doHandleSnapProblems)
     {
-        //Info<< nl
-        //    << "Introducing baffles to block off problem cells" << nl
-        //    << "----------------------------------------------" << nl
-        //    << endl;
-        //
-        //labelList facePatch
-        //(
-        //    markFacesOnProblemCells
-        //    (
-        //        motionDict,
-        //        removeEdgeConnectedCells,
-        //        perpendicularAngle,
-        //        globalToMasterPatch
-        //    )
-        //    //markFacesOnProblemCellsGeometric(motionDict)
-        //);
-        //Info<< "Analyzed problem cells in = "
-        //    << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
-        //
-        //if (debug&meshRefinement::MESH)
-        //{
-        //    faceSet problemTopo(mesh_, "problemFacesTopo", 100);
-        //
-        //    const labelList facePatchTopo
-        //    (
-        //        markFacesOnProblemCells
-        //        (
-        //            motionDict,
-        //            removeEdgeConnectedCells,
-        //            perpendicularAngle,
-        //            globalToMasterPatch
-        //        )
-        //    );
-        //    forAll(facePatchTopo, faceI)
-        //    {
-        //        if (facePatchTopo[faceI] != -1)
-        //        {
-        //            problemTopo.insert(faceI);
-        //        }
-        //    }
-        //    problemTopo.instance() = timeName();
-        //    Pout<< "Dumping " << problemTopo.size()
-        //        << " problem faces to " << problemTopo.objectPath() << endl;
-        //    problemTopo.write();
-        //}
-        //
-        //Info<< "Introducing baffles to delete problem cells." << nl << endl;
-        //
-        //if (debug)
-        //{
-        //    runTime++;
-        //}
-        //
-        //// Create baffles with same owner and neighbour for now.
-        //createBaffles(facePatch, facePatch);
-        //
-        //if (debug)
-        //{
-        //    // Debug:test all is still synced across proc patches
-        //    checkData();
-        //}
-        //Info<< "Created baffles in = "
-        //    << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
-        //
-        //printMeshInfo(debug, "After introducing baffles");
-        //
-        //if (debug&meshRefinement::MESH)
-        //{
-        //    Pout<< "Writing extra baffled mesh to time "
-        //        << timeName() << endl;
-        //    write(debug, runTime.path()/"extraBaffles");
-        //    Pout<< "Dumped debug data in = "
-        //        << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
-        //}
-
         handleSnapProblems
         (
+            snapParams,
+            useTopologicalSnapDetection,
             removeEdgeConnectedCells,
             perpendicularAngle,
             motionDict,
@@ -2192,6 +2118,8 @@ void Foam::meshRefinement::baffleAndSplitMesh
             // and delete them
             handleSnapProblems
             (
+                snapParams,
+                useTopologicalSnapDetection,
                 removeEdgeConnectedCells,
                 perpendicularAngle,
                 motionDict,

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

@@ -36,6 +36,10 @@ License
 #include "polyMeshGeometry.H"
 #include "IOmanip.H"
 #include "unitConversion.H"
+#include "autoSnapDriver.H"
+
+#include "snapParameters.H"
+#include "motionSmoother.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -968,164 +972,241 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
 }
 
 
-//// Mark faces to be baffled to prevent snapping problems. Does
+// Mark faces to be baffled to prevent snapping problems. Does
-//// test to find nearest surface and checks which faces would get squashed.
+// test to find nearest surface and checks which faces would get squashed.
-//Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
+Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
-//(
+(
-//    const dictionary& motionDict
+    const snapParameters& snapParams,
-//) const
+    const dictionary& motionDict
-//{
+) const
-//    // Construct addressing engine.
+{
-//    autoPtr<indirectPrimitivePatch> ppPtr
+    pointField oldPoints(mesh_.points());
-//    (
+
-//        meshRefinement::makePatch
+    // Repeat (most of) autoSnapDriver::doSnap
-//        (
+    {
-//            mesh_,
+        labelList adaptPatchIDs(meshedPatches());
-//            meshedPatches()
+
-//        )
+        // Construct addressing engine.
-//    );
+        autoPtr<indirectPrimitivePatch> ppPtr
-//    const indirectPrimitivePatch& pp = ppPtr();
+        (
-//    const pointField& localPoints = pp.localPoints();
+            meshRefinement::makePatch
-//    const labelList& meshPoints = pp.meshPoints();
+            (
-//
+                mesh_,
-//    // Find nearest (non-baffle) surface
+                adaptPatchIDs
-//    pointField newPoints(mesh_.points());
+            )
-//    {
+        );
-//        List<pointIndexHit> hitInfo;
+        indirectPrimitivePatch& pp = ppPtr();
-//        labelList hitSurface;
+
-//        surfaces_.findNearest
+        // Distance to attract to nearest feature on surface
-//        (
+        const scalarField snapDist
-//            surfaces_.getUnnamedSurfaces(),
+        (
-//            localPoints,
+            autoSnapDriver::calcSnapDistance(mesh_, snapParams, pp)
-//            scalarField(localPoints.size(), sqr(GREAT)),// sqr of attraction
+        );
-//            hitSurface,
+
-//            hitInfo
+
-//        );
+        // Construct iterative mesh mover.
-//
+        Info<< "Constructing mesh displacer ..." << endl;
-//        forAll(hitInfo, i)
+        Info<< "Using mesh parameters " << motionDict << nl << endl;
-//        {
+
-//            if (hitInfo[i].hit())
+        const pointMesh& pMesh = pointMesh::New(mesh_);
-//            {
+
-//                //label pointI = meshPoints[i];
+        motionSmoother meshMover
-//                //Pout<< "   " << pointI << " moved from "
+        (
-//                //    << mesh_.points()[pointI] << " by "
+            mesh_,
-//                //    << mag(hitInfo[i].hitPoint()-mesh_.points()[pointI])
+            pp,
-//                //    << endl;
+            adaptPatchIDs,
-//                newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
+            meshRefinement::makeDisplacementField(pMesh, adaptPatchIDs)(),
-//            }
+            motionDict
-//        }
+        );
-//    }
+
-//
+
-//    // Per face (internal or coupled!) the patch that the
+        // Check initial mesh
-//    // baffle should get (or -1).
+        Info<< "Checking initial mesh ..." << endl;
-//    labelList facePatch(mesh_.nFaces(), -1);
+        labelHashSet wrongFaces(mesh_.nFaces()/100);
-//    // Count of baffled faces
+        motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
-//    label nBaffleFaces = 0;
+        const label nInitErrors = returnReduce
-//
+        (
-//    {
+            wrongFaces.size(),
-//        pointField oldPoints(mesh_.points());
+            sumOp<label>()
-//        mesh_.movePoints(newPoints);
+        );
-//        faceSet wrongFaces(mesh_, "wrongFaces", 100);
+
-//        {
+        Info<< "Detected " << nInitErrors << " illegal faces"
-//            //motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
+            << " (concave, zero area or negative cell pyramid volume)"
-//
+            << endl;
-//            // Just check the errors from squashing
+
-//            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
-//
+        Info<< "Checked initial mesh in = "
-//            const labelList allFaces(identity(mesh_.nFaces()));
+            << mesh_.time().cpuTimeIncrement() << " s\n" << nl << endl;
-//            label nWrongFaces = 0;
+
-//
+        // Pre-smooth patch vertices (so before determining nearest)
-//            scalar minArea(readScalar(motionDict.lookup("minArea")));
+        autoSnapDriver::preSmoothPatch
-//            if (minArea > -SMALL)
+        (
-//            {
+            *this,
-//                polyMeshGeometry::checkFaceArea
+            snapParams,
-//                (
+            nInitErrors,
-//                    false,
+            List<labelPair>(0), //baffles
-//                    minArea,
+            meshMover
-//                    mesh_,
+        );
-//                    mesh_.faceAreas(),
+
-//                    allFaces,
+        const vectorField disp
-//                    &wrongFaces
+        (
-//                );
+            autoSnapDriver::calcNearestSurface
-//
+            (
-//                label nNewWrongFaces = returnReduce
+                *this,
-//                (
+                snapDist,   // attraction
-//                    wrongFaces.size(),
+                pp
-//                    sumOp<label>()
+            )
-//                );
+        );
-//
+
-//                Info<< "    faces with area < "
+        const labelList& meshPoints = pp.meshPoints();
-//                    << setw(5) << minArea
+
-//                    << " m^2                            : "
+        pointField newPoints(mesh_.points());
-//                    << nNewWrongFaces-nWrongFaces << endl;
+        forAll(meshPoints, i)
-//
+        {
-//                nWrongFaces = nNewWrongFaces;
+            newPoints[meshPoints[i]] += disp[i];
-//            }
+        }
-//
+        mesh_.movePoints(newPoints);
-////            scalar minDet(readScalar(motionDict.lookup("minDeterminant")));
+    }
-//            scalar minDet = 0.01;
+
-//            if (minDet > -1)
+
-//            {
+    // Per face (internal or coupled!) the patch that the
-//                polyMeshGeometry::checkCellDeterminant
+    // baffle should get (or -1).
-//                (
+    labelList facePatch(mesh_.nFaces(), -1);
-//                    false,
+    // Count of baffled faces
-//                    minDet,
+    label nBaffleFaces = 0;
-//                    mesh_,
+
-//                    mesh_.faceAreas(),
+    {
-//                    allFaces,
+        faceSet wrongFaces(mesh_, "wrongFaces", 100);
-//                    polyMeshGeometry::affectedCells(mesh_, allFaces),
+        {
-//                    &wrongFaces
+            //motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
-//                );
+
-//
+            // Just check the errors from squashing
-//                label nNewWrongFaces = returnReduce
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-//                (
+
-//                    wrongFaces.size(),
+            const labelList allFaces(identity(mesh_.nFaces()));
-//                    sumOp<label>()
+            label nWrongFaces = 0;
-//                );
+
-//
+            //const scalar minV(readScalar(motionDict.lookup("minVol", true)));
-//                Info<< "    faces on cells with determinant < "
+            //if (minV > -GREAT)
-//                    << setw(5) << minDet << "                : "
+            //{
-//                    << nNewWrongFaces-nWrongFaces << endl;
+            //    polyMeshGeometry::checkFacePyramids
-//
+            //    (
-//                nWrongFaces = nNewWrongFaces;
+            //        false,
-//            }
+            //        minV,
-//        }
+            //        mesh_,
-//
+            //        mesh_.cellCentres(),
-//
+            //        mesh_.points(),
-//        forAllConstIter(faceSet, wrongFaces, iter)
+            //        allFaces,
-//        {
+            //        List<labelPair>(0),
-//            label patchI = mesh_.boundaryMesh().whichPatch(iter.key());
+            //        &wrongFaces
-//
+            //    );
-//            if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
+            //
-//            {
+            //    label nNewWrongFaces = returnReduce
-//                facePatch[iter.key()] = getBafflePatch(facePatch, iter.key());
+            //    (
-//                nBaffleFaces++;
+            //        wrongFaces.size(),
-//
+            //        sumOp<label>()
-//                //Pout<< "    " << iter.key()
+            //    );
-//                //    //<< " on patch " << mesh_.boundaryMesh()[patchI].name()
+            //
-//                //    << " is destined for patch " << facePatch[iter.key()]
+            //    Info<< "    faces with pyramid volume < "
-//                //    << endl;
+            //        << setw(5) << minV
-//            }
+            //        << " m^3                  : "
-//        }
+            //        << nNewWrongFaces-nWrongFaces << endl;
-//        // Restore points.
+            //
-//        mesh_.movePoints(oldPoints);
+            //    nWrongFaces = nNewWrongFaces;
-//    }
+            //}
-//
+
-//
+            scalar minArea(readScalar(motionDict.lookup("minArea")));
-//    Info<< "markFacesOnProblemCellsGeometric : marked "
+            if (minArea > -SMALL)
-//        << returnReduce(nBaffleFaces, sumOp<label>())
+            {
-//        << " additional internal and coupled faces"
+                polyMeshGeometry::checkFaceArea
-//        << " to be converted into baffles." << endl;
+                (
-//
+                    false,
-//    syncTools::syncFaceList
+                    minArea,
-//    (
+                    mesh_,
-//        mesh_,
+                    mesh_.faceAreas(),
-//        facePatch,
+                    allFaces,
-//        maxEqOp<label>()
+                    &wrongFaces
-//    );
+                );
-//
+
-//    return facePatch;
+                label nNewWrongFaces = returnReduce
-//}
+                (
+                    wrongFaces.size(),
+                    sumOp<label>()
+                );
+
+                Info<< "    faces with area < "
+                    << setw(5) << minArea
+                    << " m^2                            : "
+                    << nNewWrongFaces-nWrongFaces << endl;
+
+                nWrongFaces = nNewWrongFaces;
+            }
+
+            scalar minDet(readScalar(motionDict.lookup("minDeterminant")));
+            if (minDet > -1)
+            {
+                polyMeshGeometry::checkCellDeterminant
+                (
+                    false,
+                    minDet,
+                    mesh_,
+                    mesh_.faceAreas(),
+                    allFaces,
+                    polyMeshGeometry::affectedCells(mesh_, allFaces),
+                    &wrongFaces
+                );
+
+                label nNewWrongFaces = returnReduce
+                (
+                    wrongFaces.size(),
+                    sumOp<label>()
+                );
+
+                Info<< "    faces on cells with determinant < "
+                    << setw(5) << minDet << "                : "
+                    << nNewWrongFaces-nWrongFaces << endl;
+
+                nWrongFaces = nNewWrongFaces;
+            }
+        }
+
+
+        forAllConstIter(faceSet, wrongFaces, iter)
+        {
+            label patchI = mesh_.boundaryMesh().whichPatch(iter.key());
+
+            if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
+            {
+                facePatch[iter.key()] = getBafflePatch(facePatch, iter.key());
+                nBaffleFaces++;
+
+                //Pout<< "    " << iter.key()
+                //    //<< " on patch " << mesh_.boundaryMesh()[patchI].name()
+                //    << " is destined for patch " << facePatch[iter.key()]
+                //    << endl;
+            }
+        }
+    }
+
+
+    // Restore points.
+    mesh_.movePoints(oldPoints);
+
+
+    Info<< "markFacesOnProblemCellsGeometric : marked "
+        << returnReduce(nBaffleFaces, sumOp<label>())
+        << " additional internal and coupled faces"
+        << " to be converted into baffles." << endl;
+
+    syncTools::syncFaceList
+    (
+        mesh_,
+        facePatch,
+        maxEqOp<label>()
+    );
+
+    return facePatch;
+}
 
 
 // ************************************************************************* //

src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C

@@ -77,7 +77,8 @@ const Foam::NamedEnum<Foam::refinementSurfaces::faceZoneType, 3>
 Foam::refinementSurfaces::refinementSurfaces
 (
     const searchableSurfaces& allGeometry,
-    const dictionary& surfacesDict
+    const dictionary& surfacesDict,
+    const label gapLevelIncrement
 )
 :
     allGeometry_(allGeometry),
@@ -143,7 +144,7 @@ Foam::refinementSurfaces::refinementSurfaces
             globalLevelIncr[surfI] = dict.lookupOrDefault
             (
                 "gapLevelIncrement",
-                0
+                gapLevelIncrement
             );
 
             if
@@ -274,7 +275,7 @@ Foam::refinementSurfaces::refinementSurfaces
                         label levelIncr = regionDict.lookupOrDefault
                         (
                             "gapLevelIncrement",
-                            0
+                            gapLevelIncrement
                         );
                         regionLevelIncr[surfI].insert(regionI, levelIncr);
 

src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H

@@ -149,7 +149,8 @@ public:
         refinementSurfaces
         (
             const searchableSurfaces& allGeometry,
-            const dictionary&
+            const dictionary&,
+            const label gapLevelIncrement
         );
 
         //- Construct from components

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;

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"

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

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();
 };
 
 

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_);
+    setMaster();
-    const scalar Cmu75 = pow(Cmu_, 0.75);
 
-    const scalarField& y = turbulence.y()[patchI];
+    if (patch().index() == master_)
-
-    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)
     {
-        label faceCellI = patch().faceCells()[faceI];
+        createAveragingWeights();
-
+        calculateTurbulenceFields(turbulence, G(true), epsilon(true));
-        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]);
     }
 
-    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);
 }

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
+            //- Update the coefficients associated with the patch field
-            virtual void evaluate(const Pstream::commsTypes);
+            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

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");
+        db().lookupObject<turbulenceModel>(turbulenceModel::typeName);
-    const scalarField& y = turbulence.y()[patch().index()];
 
-    const scalar Cmu25 = pow025(Cmu_);
+    setMaster();
 
-    volScalarField& G = const_cast<volScalarField&>
+    if (patch().index() == master_)
-        (
-            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)
     {
-        label faceCellI = patch().faceCells()[faceI];
+        createAveragingWeights();
-
+        calculateTurbulenceFields(turbulence, G(true), omega(true));
-        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]);
     }
 
-    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);
 }

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
 :
-    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
 

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"

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

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:
             );
         }
 
-
+    //- Destructor
-    // Member functions
+    virtual ~epsilonLowReWallFunctionFvPatchScalarField()
-
+    {}
-        // Evaluation functions
-
-            //- Update the coefficients associated with the patch field
-            virtual void updateCoeffs();
 };
 
 

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  * * * * * * * * * * * * * * //
+void epsilonWallFunctionFvPatchScalarField::setMaster()
-
-epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedInternalValueFvPatchField<scalar>(p, iF),
-    Cmu_(0.09),
-    kappa_(0.41),
-    E_(9.8)
 {
-    checkType();
+    if (master_ != -1)
-}
-
-
-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())
     {
         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] =
+        epsilon[cellI] += w*Cmu75*pow(k[cellI], 1.5)/(kappa_*y[faceI]);
-            (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
 }
 
 
-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);
 }

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
+            //- Update the coefficients associated with the patch field
-            virtual void evaluate(const Pstream::commsTypes);
+            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

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  * * * * * * * * * * * * * * //
+void omegaWallFunctionFvPatchScalarField::setMaster()
-
-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_))
 {
-    checkType();
+    if (master_ != -1)
-}
-
-
-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())
     {
         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] =
+        G[cellI] =
-            (nutw[faceI] + nuw[faceI])
+            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);
 }

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
 :
-    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
 

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;
+            }
         }
     }
 }

tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/system/controlDict

@@ -50,6 +50,7 @@ adjustTimeStep  yes;
 maxCo           0.25;
 
 maxDeltaT       1;
+maxAlphaCo      1;
 
 
 // ************************************************************************* //

tutorials/multiphase/compressibleInterFoam/laminar/depthCharge3D/system/controlDict

@@ -50,6 +50,7 @@ adjustTimeStep  yes;
 maxCo           0.25;
 
 maxDeltaT       1;
+maxAlphaCo      1;
 
 
 // ************************************************************************* //

tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phase/0.org/U

@@ -38,7 +38,7 @@ boundaryField
     }
     atmosphere
     {
-        type            fluxCorrectedVelocity;
+        type            pressureInletOutletVelocity;
         value           uniform (0 0 0);
     }
     defaultFaces

tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phaseFine/0.org/U

@@ -38,7 +38,7 @@ boundaryField
     }
     atmosphere
     {
-        type            fluxCorrectedVelocity;
+        type            pressureInletOutletVelocity;
         value           uniform (0 0 0);
     }
     defaultFaces