Merge commit 'OpenCFD/master' into olesenm

applications/solvers/incompressible/boundaryFoam/boundaryFoam.C

@@ -26,8 +26,8 @@ Application
     boundaryFoam
 
 Description
-    Steady-state solver for 1D turbulent flow, typically to generate boundary 
+    Steady-state solver for 1D turbulent flow, typically to generate boundary
-    layer conditions at an inlet, for use in a simulation. 
+    layer conditions at an inlet, for use in a simulation.
 
     Boundary layer code to calculate the U, k and epsilon distributions.
     Used to create inlet boundary conditions for experimental comparisons
@@ -82,11 +82,14 @@ int main(int argc, char *argv[])
         U += (Ubar - UbarStar);
         gradP += (Ubar - UbarStar)/(1.0/UEqn.A())().weightedAverage(mesh.V());
 
+        label id = y.size() - 1;
+
         scalar wallShearStress =
-            flowDirection & turbulence->R()()[0] & wallNormal;
+            flowDirection & turbulence->R()()[id] & wallNormal;
 
         scalar yplusWall
-            = ::sqrt(mag(wallShearStress))*y[0]/laminarTransport.nu()()[0];
+//            = ::sqrt(mag(wallShearStress))*y[id]/laminarTransport.nu()()[id];
+            = ::sqrt(mag(wallShearStress))*y[id]/turbulence->nuEff()()[id];
 
         Info<< "Uncorrected Ubar = " << (flowDirection & UbarStar.value())<< tab
             << "pressure gradient = " << (flowDirection & gradP.value()) << tab

applications/solvers/incompressible/simpleFoam/Make/options

@@ -7,5 +7,4 @@ EXE_LIBS = \
     -lincompressibleRASModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
-    -lmeshTools \
+    -lmeshTools
-    /* $(LIB_WM_OPTIONS_DIR)/libfbsdmalloc.o */

applications/utilities/postProcessing/patch/patchIntegrate/patchIntegrate.C

@@ -75,7 +75,7 @@ int main(int argc, char *argv[])
             }
 
             // Give patch area
-            Info<< "    Patch area = " << sum(mesh.Sf()) << endl;
+            Info<< "    Patch area = " << sum(mesh.Sf().boundaryField()[patchi]) << endl;
 
             if (fieldHeader.headerClassName() == "volScalarField")
             {

applications/utilities/preProcessing/applyWallFunctionBounaryConditions/Make/files

@@ -0,0 +1,3 @@
+applyWallFunctionBounaryConditions.C
+
+EXE = $(FOAM_APPBIN)/applyWallFunctionBounaryConditions

applications/utilities/preProcessing/applyWallFunctionBounaryConditions/Make/options

@@ -0,0 +1,6 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume
+

applications/utilities/preProcessing/applyWallFunctionBounaryConditions/applyWallFunctionBounaryConditions.C

@@ -0,0 +1,279 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    applyWallFunctionBounaryConditions
+
+Description
+    Updates OpenFOAM RAS cases to use the new wall function framework
+    Attempts to determine whether case is compressible or incompressible, or
+    can be supplied with -compressible command line argument
+
+\*---------------------------------------------------------------------------*/
+
+
+#include "argList.H"
+#include "fvMesh.H"
+#include "Time.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+
+#include "wallPolyPatch.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+bool caseIsCompressible(const fvMesh& mesh)
+{
+    // Attempt flux field
+    IOobject phiHeader
+    (
+        "phi",
+        mesh.time().timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    );
+
+    if (phiHeader.headerOk())
+    {
+        surfaceScalarField phi(phiHeader, mesh);
+        if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
+        {
+            return true;
+        }
+    }
+
+    // Attempt density field
+    IOobject rhoHeader
+    (
+        "rho",
+        mesh.time().timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    );
+
+    if (rhoHeader.headerOk())
+    {
+        volScalarField rho(rhoHeader, mesh);
+        if (rho.dimensions() == dimDensity)
+        {
+            return true;
+        }
+    }
+
+    // Attempt pressure field
+    IOobject pHeader
+    (
+        "p",
+        mesh.time().timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    );
+
+    if (pHeader.headerOk())
+    {
+        volScalarField p(pHeader, mesh);
+        if (p.dimensions() == dimMass/sqr(dimTime)/dimLength)
+        {
+            return true;
+        }
+    }
+
+    // Attempt hydrostatic pressure field
+    IOobject pdHeader
+    (
+        "pd",
+        mesh.time().timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    );
+
+    if (pdHeader.headerOk())
+    {
+        volScalarField pd(pdHeader, mesh);
+        if (pd.dimensions() == dimMass/sqr(dimTime)/dimLength)
+        {
+            return true;
+        }
+    }
+
+    // If none of the above are true, assume that the case is incompressible
+    return false;
+}
+
+
+void createVolScalarField
+(
+    const fvMesh& mesh,
+    const word& fieldName,
+    const dimensionSet& dims
+)
+{
+    IOobject fieldHeader
+    (
+        fieldName,
+        mesh.time().timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    );
+
+    if (!fieldHeader.headerOk())
+    {
+        Info<< "Creating field " << fieldName << nl << endl;
+
+        volScalarField field
+        (
+            IOobject
+            (
+                fieldName,
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensionedScalar("zero", dims, 0.0)
+        );
+
+        field.write();
+    }
+}
+
+
+void replaceBoundaryType
+(
+    const fvMesh& mesh,
+    const word& fieldName,
+    const word& boundaryType,
+    const string& boundaryValue
+)
+{
+    IOobject header
+    (
+        fieldName,
+        mesh.time().timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    );
+
+    if (!header.headerOk())
+    {
+        return;
+    }
+
+    Info<< "Updating boundary types for field " << header.name() << endl;
+
+    const word oldTypeName = IOdictionary::typeName;
+    const_cast<word&>(IOdictionary::typeName) = word::null;
+
+    IOdictionary dict(header);
+
+    const_cast<word&>(IOdictionary::typeName) = oldTypeName;
+    const_cast<word&>(dict.type()) = dict.headerClassName();
+
+    // Make a backup of the old field
+    word backupName(dict.name() + ".old");
+    Info<< "    copying " << dict.name() << " to "
+        << backupName << endl;
+    IOdictionary dictOld = dict;
+    dictOld.rename(backupName);
+    dictOld.regIOobject::write();
+
+    // Loop through boundary patches and update
+    const polyBoundaryMesh& bMesh = mesh.boundaryMesh();
+    dictionary& boundaryDict = dict.subDict("boundaryField");
+    forAll(bMesh, patchI)
+    {
+        if (isType<wallPolyPatch>(bMesh[patchI]))
+        {
+            word patchName = bMesh[patchI].name();
+            dictionary& oldPatch = boundaryDict.subDict(patchName);
+
+            dictionary newPatch(dictionary::null);
+            newPatch.add("type", boundaryType);
+            newPatch.add("value", ("uniform " + boundaryValue).c_str());
+
+            oldPatch = newPatch;
+        }
+    }
+
+    Info<< "    writing updated " << dict.name() << nl << endl;
+    dict.regIOobject::write();
+}
+
+
+int main(int argc, char *argv[])
+{
+
+#   include "addTimeOptions.H"
+    argList::validOptions.insert("compressible", "");
+
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createMesh.H"
+
+    bool compressible = args.options().found("compressible");
+
+    Info<< "Updating turbulence fields to operate using new run time "
+        << "selectable" << nl << "wall functions"
+        << nl << endl;
+
+    if (compressible || caseIsCompressible(mesh))
+    {
+        Info<< "Case treated as compressible" << nl << endl;
+        createVolScalarField
+        (
+            mesh,
+            "mut",
+            dimArea/dimTime*dimDensity
+        );
+        replaceBoundaryType(mesh, "mut", "mutWallFunction", "0");
+    }
+    else
+    {
+        Info<< "Case treated as incompressible" << nl << endl;
+        createVolScalarField(mesh, "nut", dimArea/dimTime);
+        replaceBoundaryType(mesh, "nut", "nutWallFunction", "0");
+    }
+
+    replaceBoundaryType(mesh, "epsilon", "epsilonWallFunction", "0");
+    replaceBoundaryType(mesh, "omega", "omegaWallFunction", "0");
+    replaceBoundaryType(mesh, "k", "kQRWallFunction", "0");
+    replaceBoundaryType(mesh, "q", "kQRWallFunction", "0");
+    replaceBoundaryType(mesh, "R", "kQRWallFunction", "(0 0 0 0 0 0)");
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

src/OSspecific/Unix/regularExpression.H

@@ -0,0 +1,121 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::regularExpression
+
+Description
+    Wrapper around regular expressions.
+
+SourceFiles
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef regularExpression_H
+#define regularExpression_H
+
+#include <sys/types.h>
+#include <regex.h>
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class regularExpression Declaration
+\*---------------------------------------------------------------------------*/
+
+class regularExpression
+{
+    // Private data
+
+        //- Precompiled regular expression
+        regex_t* preg_;
+
+
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        regularExpression(const regularExpression&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const regularExpression&);
+
+public:
+
+
+    // Constructors
+
+        //- Construct from string
+        inline regularExpression(const string& s)
+        {
+            preg_ = new regex_t;
+
+            if (regcomp(preg_, s.c_str(), REG_EXTENDED|REG_NOSUB) != 0)
+            {
+                FatalErrorIn
+                (
+                    "regularExpression::regularExpression(const char*)"
+                )   << "Failed to compile regular expression " << s
+                    << exit(FatalError);
+            }
+        }
+
+
+    // Destructor
+
+        //- Construct from string
+        inline ~regularExpression()
+        {
+            if (preg_)
+            {
+                regfree(preg_);
+                delete preg_;
+            }
+        }
+
+
+    // Member functions
+
+        //- Matches?
+        inline bool matches(const string& s)
+        {
+            size_t nmatch = 0;
+            regmatch_t *pmatch = NULL;
+
+            return regexec(preg_, s.c_str(), nmatch, pmatch, 0) == 0;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/containers/Lists/DynamicList/DynamicList.H

@@ -159,15 +159,9 @@ public:
         //- Assignment of all entries to the given value
         inline void operator=(const T&);
 
-        //- Assignment from List<T>
+        //- Assignment from List<T>. Also handles assignment from DynamicList.
         inline void operator=(const List<T>&);
 
-        //- Assignment from DynamicList<T>
-        inline void operator=
-        (
-            const DynamicList<T, SizeInc, SizeMult, SizeDiv>&
-        );
-
 
     // IOstream operators
 

src/OpenFOAM/containers/Lists/DynamicList/DynamicListI.H

@@ -85,8 +85,10 @@ inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::setSize
     }
     else
     {
+        label nextFree = List<T>::size();
         allocSize_ = s;
         List<T>::setSize(allocSize_);
+        List<T>::size() = nextFree;
     }
 }
 
@@ -104,8 +106,10 @@ inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::setSize
     }
     else
     {
+        label nextFree = List<T>::size();
         allocSize_ = s;
         List<T>::setSize(allocSize_, t);
+        List<T>::size() = nextFree;
     }
 }
 
@@ -165,7 +169,7 @@ Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::transfer
 )
 {
     allocSize_ = l.allocSize();
-    List<T>::transfer(l);       // take over storage
+    List<T>::transfer(l);       // take over storage. Null l.
 }
 
 
@@ -187,9 +191,9 @@ inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::append(const T& e)
         List<T>::setSize(allocSize_);
     }
 
-    this->operator[](nextFree - 1) = e;
-
     List<T>::size() = nextFree;
+
+    this->operator[](nextFree - 1) = e;
 }
 
 
@@ -204,7 +208,13 @@ inline T Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::remove()
         )   << "List is empty" << abort(FatalError);
     }
 
-    return List<T>::operator[](--List<T>::size());
+    label nextFree = List<T>::size()-1;
+
+    const T& val = List<T>::operator[](nextFree);
+
+    List<T>::size() = nextFree;
+
+    return val;
 }
 
 
@@ -258,18 +268,4 @@ inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::operator=
 }
 
 
-template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
-inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::operator=
-(
-    const DynamicList<T, SizeInc, SizeMult, SizeDiv>& l
-)
-{
-    List<T>::operator=(l);
-    // allocSize_ = l.allocSize();  // wrong
-    allocSize_ = List<T>::size();
-    // ^^^^ with this change, we could just use
-    //      DynamicList::operator=(const List<T>&) instead
-}
-
-
 // ************************************************************************* //

src/OpenFOAM/meshes/meshShapes/face/face.C

@@ -526,7 +526,7 @@ Foam::point Foam::face::centre(const pointField& meshPoints) const
 }
 
 
-Foam::vector Foam::face::normal(const pointField& meshPoints) const
+Foam::vector Foam::face::normal(const pointField& p) const
 {
     // Calculate the normal by summing the face triangle normals.
     // Changed to deal with small concavity by using a central decomposition
@@ -539,38 +539,43 @@ Foam::vector Foam::face::normal(const pointField& meshPoints) const
     {
         return triPointRef
         (
-            meshPoints[operator[](0)],
+            p[operator[](0)],
-            meshPoints[operator[](1)],
+            p[operator[](1)],
-            meshPoints[operator[](2)]
+            p[operator[](2)]
         ).normal();
     }
 
-    vector n = vector::zero;
-
-    point centrePoint = Foam::average(points(meshPoints));
-
     label nPoints = size();
 
-    point nextPoint = centrePoint;
-
     register label pI;
 
+    point centrePoint = vector::zero;
+    for (pI = 0; pI < nPoints; pI++)
+    {
+        centrePoint += p[operator[](pI)];
+    }
+    centrePoint /= nPoints;
+
+    vector n = vector::zero;
+
+    point nextPoint = centrePoint;
+
     for (pI = 0; pI < nPoints; pI++)
     {
         if (pI < nPoints - 1)
         {
-            nextPoint = meshPoints[operator[](pI + 1)];
+            nextPoint = p[operator[](pI + 1)];
         }
         else
         {
-            nextPoint = meshPoints[operator[](0)];
+            nextPoint = p[operator[](0)];
         }
 
         // Note: for best accuracy, centre point always comes last
         //
         n += triPointRef
         (
-            meshPoints[operator[](pI)],
+            p[operator[](pI)],
             nextPoint,
             centrePoint
         ).normal();

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C

@@ -241,7 +241,7 @@ void Foam::mapDistribute::distribute
                     Pstream::nonBlocking,
                     domain,
                     reinterpret_cast<const char*>(subField.begin()),
-                    subField.size()
+                    subField.size()*sizeof(T)
                 );
             }
         }
@@ -262,7 +262,7 @@ void Foam::mapDistribute::distribute
                     Pstream::nonBlocking,
                     domain,
                     reinterpret_cast<char*>(recvFields[domain].begin()),
-                    recvFields[domain].size()
+                    recvFields[domain].size()*sizeof(T)
                 );
             }
         }

src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C

@@ -66,7 +66,6 @@ primitiveMesh::primitiveMesh()
     ppPtr_(NULL),
     cpPtr_(NULL),
 
-    allocSize_(0),
     labels_(0),
 
     cellCentresPtr_(NULL),
@@ -109,7 +108,6 @@ primitiveMesh::primitiveMesh
     ppPtr_(NULL),
     cpPtr_(NULL),
 
-    allocSize_(0),
     labels_(0),
 
     cellCentresPtr_(NULL),

src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H

@@ -54,6 +54,7 @@ SourceFiles
 #ifndef primitiveMesh_H
 #define primitiveMesh_H
 
+#include "DynamicList.H"
 #include "edgeList.H"
 #include "pointField.H"
 #include "SubField.H"
@@ -157,10 +158,8 @@ class primitiveMesh
 
         // On-the-fly edge addresing storage
 
-            //- Temporary storage for addressing. allocSize is the real size
+            //- Temporary storage for addressing.
-            //  of the labelList.
+            mutable DynamicList<label> labels_;
-            mutable label allocSize_;
-            mutable labelList labels_;
 
             //- Temporary storage for addressing
             mutable labelHashSet labelSet_;
@@ -705,31 +704,80 @@ public:
 
             // On-the-fly addressing calculation. These functions return either
             // a reference to the full addressing (if already calculated) or
-            // a reference to member data labels_ so be careful when not storing
+            // a reference to the supplied storage. The one-argument ones
+            // use member DynamicList labels_ so be careful when not storing
             // result.
 
-            //- cellCells using cells
+            //- cellCells using cells.
+            const labelList& cellCells
+            (
+                const label cellI,
+                DynamicList<label>&
+            ) const;
+
             const labelList& cellCells(const label cellI) const;
 
             //- cellPoints using cells
+            const labelList& cellPoints
+            (
+                const label cellI,
+                DynamicList<label>&
+            ) const;
+
             const labelList& cellPoints(const label cellI) const;
 
             //- pointCells using pointFaces
+            const labelList& pointCells
+            (
+                const label pointI,
+                DynamicList<label>&
+            ) const;
+
             const labelList& pointCells(const label pointI) const;
 
             //- pointPoints using edges, pointEdges
+            const labelList& pointPoints
+            (
+                const label pointI,
+                DynamicList<label>&
+            ) const;
+
             const labelList& pointPoints(const label pointI) const;
 
             //- faceEdges using pointFaces, edges, pointEdges
+            const labelList& faceEdges
+            (
+                const label faceI,
+                DynamicList<label>&
+            ) const;
+
             const labelList& faceEdges(const label faceI) const;
 
             //- edgeFaces using pointFaces, edges, pointEdges
+            const labelList& edgeFaces
+            (
+                const label edgeI,
+                DynamicList<label>&
+            ) const;
+
             const labelList& edgeFaces(const label edgeI) const;
 
             //- edgeCells using pointFaces, edges, pointEdges
+            const labelList& edgeCells
+            (
+                const label edgeI,
+                DynamicList<label>&
+            ) const;
+
             const labelList& edgeCells(const label edgeI) const;
 
             //- cellEdges using cells, pointFaces, edges, pointEdges
+            const labelList& cellEdges
+            (
+                const label cellI,
+                DynamicList<label>&
+            ) const;
+
             const labelList& cellEdges(const label cellI) const;
 
 

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCells.C

@@ -41,6 +41,14 @@ void primitiveMesh::calcCellCells() const
     {
         Pout<< "primitiveMesh::calcCellCells() : calculating cellCells"
             << endl;
+
+        if (debug == -1)
+        {
+            // For checking calls:abort so we can quickly hunt down
+            // origin of call
+            FatalErrorIn("primitiveMesh::calcCellCells()")
+                << abort(FatalError);
+        }
     }
 
     // It is an error to attempt to recalculate cellCells
@@ -105,7 +113,11 @@ const labelListList& primitiveMesh::cellCells() const
 }
 
 
-const labelList& primitiveMesh::cellCells(const label cellI) const
+const labelList& primitiveMesh::cellCells
+(
+    const label cellI,
+    DynamicList<label>& storage
+) const
 {
     if (hasCellCells())
     {
@@ -117,16 +129,7 @@ const labelList& primitiveMesh::cellCells(const label cellI) const
         const labelList& nei = faceNeighbour();
         const cell& cFaces = cells()[cellI];
 
-        labels_.size() = allocSize_;
+        storage.clear();
-
-        if (cFaces.size() > allocSize_)
-        {
-            labels_.clear();
-            allocSize_ = cFaces.size();
-            labels_.setSize(allocSize_);
-        }
-
-        label n = 0;
 
         forAll(cFaces, i)
         {
@@ -136,22 +139,26 @@ const labelList& primitiveMesh::cellCells(const label cellI) const
             {
                 if (own[faceI] == cellI)
                 {
-                    labels_[n++] = nei[faceI];
+                    storage.append(nei[faceI]);
                 }
                 else
                 {
-                    labels_[n++] = own[faceI];
+                    storage.append(own[faceI]);
                 }
             }
         }
 
-        labels_.size() = n;
+        return storage;
-
-        return labels_;
     }
 }
 
 
+const labelList& primitiveMesh::cellCells(const label cellI) const
+{
+    return cellCells(cellI, labels_);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellEdges.C

@@ -40,6 +40,14 @@ void Foam::primitiveMesh::calcCellEdges() const
         Pout<< "primitiveMesh::calcCellEdges() : "
             << "calculating cellEdges"
             << endl;
+
+        if (debug == -1)
+        {
+            // For checking calls:abort so we can quickly hunt down
+            // origin of call
+            FatalErrorIn("primitiveMesh::calcCellEdges()")
+                << abort(FatalError);
+        }
     }
 
     // It is an error to attempt to recalculate cellEdges

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellPoints.C

@@ -42,6 +42,14 @@ const labelListList& primitiveMesh::cellPoints() const
         {
             Pout<< "primitiveMesh::cellPoints() : "
                 << "calculating cellPoints" << endl;
+
+            if (debug == -1)
+            {
+                // For checking calls:abort so we can quickly hunt down
+                // origin of call
+                FatalErrorIn("primitiveMesh::cellPoints()")
+                    << abort(FatalError);
+            }
         }
 
         // Invert pointCells
@@ -53,7 +61,11 @@ const labelListList& primitiveMesh::cellPoints() const
 }
 
 
-const labelList& primitiveMesh::cellPoints(const label cellI) const
+const labelList& primitiveMesh::cellPoints
+(
+    const label cellI,
+    DynamicList<label>& storage
+) const
 {
     if (hasCellPoints())
     {
@@ -76,29 +88,28 @@ const labelList& primitiveMesh::cellPoints(const label cellI) const
             }
         }
 
-        labels_.size() = allocSize_;
+        storage.clear();
-
+        if (labelSet_.size() > storage.allocSize())
-        if (labelSet_.size() > allocSize_)
         {
-            labels_.clear();
+            storage.setSize(labelSet_.size());
-            allocSize_ = labelSet_.size();
-            labels_.setSize(allocSize_);
         }
 
-        label n = 0;
-
         forAllConstIter(labelHashSet, labelSet_, iter)
         {
-            labels_[n++] = iter.key();
+            storage.append(iter.key());
         }
 
-        labels_.size() = n;
+        return storage;
-
-        return labels_;
     }
 }
 
 
+const labelList& primitiveMesh::cellPoints(const label cellI) const
+{
+    return cellPoints(cellI, labels_);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdgeCells.C

@@ -41,6 +41,14 @@ const labelListList& primitiveMesh::edgeCells() const
         if (debug)
         {
             Pout<< "primitiveMesh::edgeCells() : calculating edgeCells" << endl;
+
+            if (debug == -1)
+            {
+                // For checking calls:abort so we can quickly hunt down
+                // origin of call
+                FatalErrorIn("primitiveMesh::edgeCells()")
+                    << abort(FatalError);
+            }
         }
         // Invert cellEdges
         ecPtr_ = new labelListList(nEdges());
@@ -51,7 +59,11 @@ const labelListList& primitiveMesh::edgeCells() const
 }
 
 
-const labelList& primitiveMesh::edgeCells(const label edgeI) const
+const labelList& primitiveMesh::edgeCells
+(
+    const label edgeI,
+    DynamicList<label>& storage
+) const
 {
     if (hasEdgeCells())
     {
@@ -62,24 +74,11 @@ const labelList& primitiveMesh::edgeCells(const label edgeI) const
         const labelList& own = faceOwner();
         const labelList& nei = faceNeighbour();
 
-        // edge faces can either return labels_ or reference in edgeLabels.
+        // Construct edgeFaces
-        labelList labelsCopy;
+        DynamicList<label> eFacesStorage;
-        if (!hasEdgeFaces())
+        const labelList& eFaces = edgeFaces(edgeI, eFacesStorage);
-        {
-            labelsCopy = edgeFaces(edgeI);
-        }
 
-        const labelList& eFaces =
+        storage.clear();
-        (
-            hasEdgeFaces()
-          ? edgeFaces()[edgeI]
-          : labelsCopy
-        );
-
-        labels_.size() = allocSize_;
-
-        // labels_ should certainly be big enough for edge cells.
-        label n = 0;
 
         // Do quadratic insertion.
         forAll(eFaces, i)
@@ -89,10 +88,10 @@ const labelList& primitiveMesh::edgeCells(const label edgeI) const
             {
                 label ownCellI = own[faceI];
 
-                // Check if not already in labels_
+                // Check if not already in storage
-                for (label j = 0; j < n; j++)
+                forAll(storage, j)
                 {
-                    if (labels_[j] == ownCellI)
+                    if (storage[j] == ownCellI)
                     {
                         ownCellI = -1;
                         break;
@@ -101,7 +100,7 @@ const labelList& primitiveMesh::edgeCells(const label edgeI) const
 
                 if (ownCellI != -1)
                 {
-                    labels_[n++] = ownCellI;
+                    storage.append(ownCellI);
                 }
             }
 
@@ -109,9 +108,9 @@ const labelList& primitiveMesh::edgeCells(const label edgeI) const
             {
                 label neiCellI = nei[faceI];
 
-                for (label j = 0; j < n; j++)
+                forAll(storage, j)
                 {
-                    if (labels_[j] == neiCellI)
+                    if (storage[j] == neiCellI)
                     {
                         neiCellI = -1;
                         break;
@@ -120,18 +119,22 @@ const labelList& primitiveMesh::edgeCells(const label edgeI) const
 
                 if (neiCellI != -1)
                 {
-                    labels_[n++] = neiCellI;
+                    storage.append(neiCellI);
                 }
             }
         }
 
-        labels_.size() = n;
+        return storage;
-
-        return labels_;
     }
 }
 
 
+const labelList& primitiveMesh::edgeCells(const label edgeI) const
+{
+    return edgeCells(edgeI, labels_);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdgeFaces.C

@@ -41,6 +41,14 @@ const labelListList& primitiveMesh::edgeFaces() const
         if (debug)
         {
             Pout<< "primitiveMesh::edgeFaces() : calculating edgeFaces" << endl;
+
+            if (debug == -1)
+            {
+                // For checking calls:abort so we can quickly hunt down
+                // origin of call
+                FatalErrorIn("primitiveMesh::edgeFaces()")
+                    << abort(FatalError);
+            }
         }
 
         // Invert faceEdges
@@ -52,7 +60,11 @@ const labelListList& primitiveMesh::edgeFaces() const
 }
 
 
-const labelList& primitiveMesh::edgeFaces(const label edgeI) const
+const labelList& primitiveMesh::edgeFaces
+(
+    const label edgeI,
+    DynamicList<label>& storage
+) const
 {
     if (hasEdgeFaces())
     {
@@ -67,9 +79,8 @@ const labelList& primitiveMesh::edgeFaces(const label edgeI) const
 
         label i0 = 0;
         label i1 = 0;
-        label n = 0;
 
-        labels_.size() = allocSize_;
+        storage.clear();
 
         while (i0 < pFaces0.size() && i1 < pFaces1.size())
         {
@@ -84,26 +95,23 @@ const labelList& primitiveMesh::edgeFaces(const label edgeI) const
             else
             {
                 // Equal. Append.
-                if (n == allocSize_)
+                storage.append(pFaces0[i0]);
-                {
-                    // Have setSize copy contents so far
-                    labels_.size() = n;
-                    allocSize_ = allocSize_*2 + 1;
-                    labels_.setSize(allocSize_);
-                }
-                labels_[n++] = pFaces0[i0];
                 ++i0;
                 ++i1;
             }
         }
 
-        labels_.size() = n;
+        return storage;
-
-        return labels_;
     }
 }
 
 
+const labelList& primitiveMesh::edgeFaces(const label edgeI) const
+{
+    return edgeFaces(edgeI, labels_);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdges.C

@@ -582,11 +582,15 @@ void primitiveMesh::clearOutEdges()
     deleteDemandDrivenData(pePtr_);
     deleteDemandDrivenData(fePtr_);
     labels_.clear();
-    allocSize_ = 0;
+    labelSet_.clear();
 }
 
 
-const labelList& primitiveMesh::faceEdges(const label faceI) const
+const labelList& primitiveMesh::faceEdges
+(
+    const label faceI,
+    DynamicList<label>& storage
+) const
 {
     if (hasFaceEdges())
     {
@@ -597,34 +601,40 @@ const labelList& primitiveMesh::faceEdges(const label faceI) const
         const labelListList& pointEs = pointEdges();
         const face& f = faces()[faceI];
 
-        labels_.size() = allocSize_;
+        storage.clear();
-
+        if (f.size() > storage.allocSize())
-        if (f.size() > allocSize_)
         {
-            labels_.clear();
+            storage.setSize(f.size());
-            allocSize_ = f.size();
-            labels_.setSize(allocSize_);
         }
 
-        label n = 0;
-
         forAll(f, fp)
         {
-            labels_[n++] = findFirstCommonElementFromSortedLists
+            storage.append
             (
-                pointEs[f[fp]],
+                findFirstCommonElementFromSortedLists
-                pointEs[f.nextLabel(fp)]
+                (
+                    pointEs[f[fp]],
+                    pointEs[f.nextLabel(fp)]
+                )
             );
         }
 
-        labels_.size() = n;
+        return storage;
-
-        return labels_;
     }
 }
 
 
-const labelList& primitiveMesh::cellEdges(const label cellI) const
+const labelList& primitiveMesh::faceEdges(const label faceI) const
+{
+    return faceEdges(faceI, labels_);
+}
+
+
+const labelList& primitiveMesh::cellEdges
+(
+    const label cellI,
+    DynamicList<label>& storage
+) const
 {
     if (hasCellEdges())
     {
@@ -646,29 +656,29 @@ const labelList& primitiveMesh::cellEdges(const label cellI) const
             }
         }
 
-        labels_.size() = allocSize_;
+        storage.clear();
 
-        if (labelSet_.size() > allocSize_)
+        if (labelSet_.size() > storage.allocSize())
         {
-            labels_.clear();
+            storage.setSize(labelSet_.size());
-            allocSize_ = labelSet_.size();
-            labels_.setSize(allocSize_);
         }
 
-        label n =0;
-
         forAllConstIter(labelHashSet, labelSet_, iter)
         {
-            labels_[n++] = iter.key();
+            storage.append(iter.key());
         }
 
-        labels_.size() = n;
+        return storage;
-
-        return labels_;
     }
 }
 
 
+const labelList& primitiveMesh::cellEdges(const label cellI) const
+{
+    return cellEdges(cellI, labels_);;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshPointCells.C

@@ -43,6 +43,14 @@ void primitiveMesh::calcPointCells() const
         Pout<< "primitiveMesh::calcPointCells() : "
             << "calculating pointCells"
             << endl;
+
+        if (debug == -1)
+        {
+            // For checking calls:abort so we can quickly hunt down
+            // origin of call
+            FatalErrorIn("primitiveMesh::calcPointCells()")
+                << abort(FatalError);
+        }
     }
 
     // It is an error to attempt to recalculate pointCells
@@ -114,7 +122,11 @@ const labelListList& primitiveMesh::pointCells() const
 }
 
 
-const labelList& primitiveMesh::pointCells(const label pointI) const
+const labelList& primitiveMesh::pointCells
+(
+    const label pointI,
+    DynamicList<label>& storage
+) const
 {
     if (hasPointCells())
     {
@@ -126,58 +138,48 @@ const labelList& primitiveMesh::pointCells(const label pointI) const
         const labelList& nei = faceNeighbour();
         const labelList& pFaces = pointFaces()[pointI];
 
-        labels_.size() = allocSize_;
+        storage.clear();
-
-        label n = 0;
 
         forAll(pFaces, i)
         {
             const label faceI = pFaces[i];
 
             // Append owner
-            if (n == allocSize_)
+            storage.append(own[faceI]);
-            {
-                labels_.size() = n;
-                allocSize_ = allocSize_*2 + 1;
-                labels_.setSize(allocSize_);
-            }
-            labels_[n++] = own[faceI];
 
             // Append neighbour
             if (faceI < nInternalFaces())
             {
-                if (n == allocSize_)
+                storage.append(nei[faceI]);
-                {
-                    labels_.size() = n;
-                    allocSize_ = allocSize_*2 + 1;
-                    labels_.setSize(allocSize_);
-                }
-                labels_[n++] = nei[faceI];
             }
         }
-        labels_.size() = n;
-
 
         // Filter duplicates
-        sort(labels_);
+        sort(storage);
 
-        n = 1;
+        label n = 1;
 
-        for (label i = 1; i < labels_.size(); i++)
+        for (label i = 1; i < storage.size(); i++)
         {
-            if (labels_[i] != labels_[i-1])
+            if (storage[i] != storage[i-1])
             {
-                labels_[n++] = labels_[i];
+                storage[n++] = storage[i];
             }
         }
 
-        labels_.size() = n;
+        storage.size() = n;
 
-        return labels_;
+        return storage;
     }
 }
 
 
+const labelList& primitiveMesh::pointCells(const label pointI) const
+{
+    return pointCells(pointI, labels_);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshPointPoints.C

@@ -40,6 +40,14 @@ void primitiveMesh::calcPointPoints() const
         Pout<< "primitiveMesh::calcPointPoints() : "
             << "calculating pointPoints"
             << endl;
+
+        if (debug == -1)
+        {
+            // For checking calls:abort so we can quickly hunt down
+            // origin of call
+            FatalErrorIn("primitiveMesh::calcPointPoints()")
+                << abort(FatalError);
+        }
     }
 
     // It is an error to attempt to recalculate pointPoints
@@ -97,7 +105,11 @@ const labelListList& primitiveMesh::pointPoints() const
 }
 
 
-const labelList& primitiveMesh::pointPoints(const label pointI) const
+const labelList& primitiveMesh::pointPoints
+(
+    const label pointI,
+    DynamicList<label>& storage
+) const
 {
     if (hasPointPoints())
     {
@@ -108,30 +120,29 @@ const labelList& primitiveMesh::pointPoints(const label pointI) const
         const edgeList& edges = this->edges();
         const labelList& pEdges = pointEdges()[pointI];
 
-        labels_.size() = allocSize_;
+        storage.clear();
 
-        if (pEdges.size() > allocSize_)
+        if (pEdges.size() > storage.allocSize())
         {
-            // Set size() so memory allocation behaves as normal.
+            storage.setSize(pEdges.size());
-            labels_.clear();
-            allocSize_ = pEdges.size();
-            labels_.setSize(allocSize_);
         }
 
-        label n = 0;
-
         forAll(pEdges, i)
         {
-            labels_[n++] = edges[pEdges[i]].otherVertex(pointI);
+            storage.append(edges[pEdges[i]].otherVertex(pointI));
         }
 
-        labels_.size() = n;
+        return storage;
-
-        return labels_;
     }
 }
 
 
+const labelList& primitiveMesh::pointPoints(const label pointI) const
+{
+    return pointPoints(pointI, labels_);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/primitives/Lists/stringListTemplates.C

@@ -25,9 +25,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "labelList.H"
-
+#include "regularExpression.H"
-#include <sys/types.h>
-#include <regex.h>
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -41,24 +39,12 @@ labelList findStrings(const string& regexp, const StringList& sl)
 {
     labelList matches(sl.size());
 
-    regex_t *preg = new regex_t;
+    regularExpression re(regexp);
-
-    if (regcomp(preg, regexp.c_str(), REG_EXTENDED|REG_NOSUB) != 0)
-    {
-        WarningIn("findStrings(const string& regexp, const stringList& sl)")
-            << "Failed to compile regular expression " << regexp
-            << endl;
-
-        return matches;
-    }
-
-    size_t nmatch = 0;
-    regmatch_t *pmatch = NULL;
 
     label matchi = 0;
     forAll(sl, i)
     {
-        if (regexec(preg, sl[i].c_str(), nmatch, pmatch, 0) == 0)
+        if (re.matches(sl[i]))
         {
             matches[matchi++] = i;
         }
@@ -66,9 +52,6 @@ labelList findStrings(const string& regexp, const StringList& sl)
 
     matches.setSize(matchi);
 
-    regfree(preg);
-    delete preg;
-
     return matches;
 }
 

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C

@@ -1059,6 +1059,8 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
         // Surfaces with zone information
         const wordList& faceZoneNames = surfaces.faceZoneNames();
 
+        scalarField minSnapDist(snapDist);
+
         forAll(zonedSurfaces, i)
         {
             label zoneSurfI = zonedSurfaces[i];
@@ -1075,31 +1077,33 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
                 )
             );
 
-            pointField zonePoints(zonePointIndices.size());
-            forAll(zonePointIndices, i)
-            {
-                zonePoints[i] = localPoints[zonePointIndices[i]];
-            }
-
             // Find nearest for points both on faceZone and pp.
             List<pointIndexHit> hitInfo;
             labelList hitSurface;
             surfaces.findNearest
             (
                 labelList(1, zoneSurfI),
-                zonePoints,
+                pointField(localPoints, zonePointIndices),
-                sqr(4*snapDist),
+                sqr(4*scalarField(minSnapDist, zonePointIndices)),
                 hitSurface,
                 hitInfo
             );
 
-            forAll(hitInfo, pointI)
+            forAll(hitInfo, i)
             {
-                if (hitInfo[pointI].hit())
+                label pointI = zonePointIndices[i]; 
+
+                if (hitInfo[i].hit())
                 {
                     patchDisp[pointI] =
-                        hitInfo[pointI].hitPoint()
+                        hitInfo[i].hitPoint()
                       - localPoints[pointI];
+
+                    minSnapDist[pointI] = min
+                    (
+                        minSnapDist[pointI],
+                        mag(patchDisp[pointI])
+                    );
                 }
                 else
                 {
@@ -1107,7 +1111,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
                         << "For point:" << pointI
                         << " coordinate:" << localPoints[pointI]
                         << " did not find any surface within:"
-                        << 4*snapDist[pointI]
+                        << 4*minSnapDist[pointI]
                         << " meter." << endl;
                 }
             }

src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C

@@ -514,7 +514,7 @@ void Foam::meshRefinement::markBoundaryFace
 {
     isBoundaryFace[faceI] = true;
 
-    const labelList& fEdges = mesh_.faceEdges()[faceI];
+    const labelList& fEdges = mesh_.faceEdges(faceI);
 
     forAll(fEdges, fp)
     {
@@ -623,12 +623,16 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
     // If so what is the remaining non-boundary anchor point?
     labelHashSet nonBoundaryAnchors(mesh_.nCells()/10000);
 
+    // On-the-fly addressing storage.
+    DynamicList<label> dynFEdges;
+    DynamicList<label> dynCPoints;
+
     // Count of faces marked for baffling
     label nBaffleFaces = 0;
 
     forAll(cellLevel, cellI)
     {
-        const labelList cPoints(meshCutter_.cellPoints(cellI));
+        const labelList& cPoints = mesh_.cellPoints(cellI, dynCPoints);
 
         // Get number of anchor points (pointLevel == cellLevel)
 
@@ -714,11 +718,14 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
     // Loop over all points. If a point is connected to 4 or more cells
     // with 7 anchor points on the boundary set those cell's non-boundary faces
     // to baffles
+
+    DynamicList<label> dynPCells;
+
     forAllConstIter(labelHashSet, nonBoundaryAnchors, iter)
     {
         label pointI = iter.key();
 
-        const labelList& pCells = mesh_.pointCells()[pointI];
+        const labelList& pCells = mesh_.pointCells(pointI, dynPCells);
 
         // Count number of 'hasSevenBoundaryAnchorPoints' cells.
         label n = 0;
@@ -806,7 +813,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
     {
         if (facePatch[faceI] == -1)
         {
-            const labelList& fEdges = mesh_.faceEdges()[faceI];
+            const labelList& fEdges = mesh_.faceEdges(faceI, dynFEdges);
             label nFaceBoundaryEdges = 0;
 
             forAll(fEdges, fe)
@@ -840,7 +847,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
             {
                 if (facePatch[faceI] == -1)
                 {
-                    const labelList& fEdges = mesh_.faceEdges()[faceI];
+                    const labelList& fEdges = mesh_.faceEdges(faceI, dynFEdges);
                     label nFaceBoundaryEdges = 0;
 
                     forAll(fEdges, fe)
@@ -1239,6 +1246,7 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::filterDuplicateFaces
     labelList nBafflesPerEdge(mesh_.nEdges(), 0);
 
 
+
     // Count number of boundary faces per edge
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -1255,7 +1263,7 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::filterDuplicateFaces
 
             forAll(pp, i)
             {
-                const labelList& fEdges = mesh_.faceEdges()[faceI];
+                const labelList& fEdges = mesh_.faceEdges(faceI);
 
                 forAll(fEdges, fEdgeI)
                 {
@@ -1267,19 +1275,23 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::filterDuplicateFaces
     }
 
 
+    DynamicList<label> fe0;
+    DynamicList<label> fe1;
+
+
     // Count number of duplicate boundary faces per edge
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     forAll(couples, i)
     {
-        const labelList& fEdges0 = mesh_.faceEdges()[couples[i].first()];
+        const labelList& fEdges0 = mesh_.faceEdges(couples[i].first(), fe0);
 
         forAll(fEdges0, fEdgeI)
         {
             nBafflesPerEdge[fEdges0[fEdgeI]] += 1000000;
         }
 
-        const labelList& fEdges1 = mesh_.faceEdges()[couples[i].second()];
+        const labelList& fEdges1 = mesh_.faceEdges(couples[i].second(), fe1);
 
         forAll(fEdges1, fEdgeI)
         {
@@ -1314,7 +1326,7 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::filterDuplicateFaces
          == patches.whichPatch(couple.second())
         )
         {
-            const labelList& fEdges = mesh_.faceEdges()[couples[i].first()];
+            const labelList& fEdges = mesh_.faceEdges(couples[i].first());
 
             forAll(fEdges, fEdgeI)
             {

src/dynamicMesh/motionSmoother/motionSmoother.C

@@ -343,7 +343,7 @@ void Foam::motionSmoother::getAffectedFacesAndPoints
 
         forAllConstIter(pointSet, nbrPoints, iter)
         {
-            const labelList& pCells = mesh_.pointCells()[iter.key()];  
+            const labelList& pCells = mesh_.pointCells(iter.key());
 
             forAll(pCells, pCellI)
             {

src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C

@@ -330,15 +330,14 @@ Foam::label Foam::addPatchCellLayer::addSideFace
     const label meshEdgeI,              // corresponding mesh edge
     const label layerI,                 // layer
     const label numEdgeFaces,           // number of layers for edge
+    const labelList& meshFaces,         // precalculated edgeFaces
     polyTopoChange& meshMod
 ) const
 {
     // Edge to 'inflate' from
     label inflateEdgeI = -1;
 
-    // Mesh faces using edge
+    // Check mesh faces using edge
-    const labelList& meshFaces = mesh_.edgeFaces()[meshEdgeI];
-
     forAll(meshFaces, i)
     {
         if (mesh_.isInternalFace(meshFaces[i]))
@@ -620,6 +619,9 @@ void Foam::addPatchCellLayer::setRefinement
 
     const labelList& meshPoints = pp.meshPoints();
 
+    // Some storage for edge-face-addressing.
+    DynamicList<label> ef;
+
     // Precalculate mesh edges for pp.edges.
     labelList meshEdges(calcMeshEdges(mesh_, pp));
 
@@ -777,7 +779,9 @@ void Foam::addPatchCellLayer::setRefinement
                 label meshEdgeI = meshEdges[edgeI];
 
                 // Mesh faces using edge
-                const labelList& meshFaces = mesh_.edgeFaces()[meshEdgeI];
+
+                // Mesh faces using edge
+                const labelList& meshFaces = mesh_.edgeFaces(meshEdgeI, ef);
 
                 // Check that there is only one patchface using edge.
                 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
@@ -1353,6 +1357,12 @@ void Foam::addPatchCellLayer::setRefinement
                             patchFaceI
                         );
 
+                        const labelList& meshFaces = mesh_.edgeFaces
+                        (
+                            meshEdgeI,
+                            ef
+                        );
+
                         addSideFace
                         (
                             pp,
@@ -1365,6 +1375,7 @@ void Foam::addPatchCellLayer::setRefinement
                             meshEdgeI,      // corresponding mesh edge
                             i,
                             numEdgeSideFaces,
+                            meshFaces,
                             meshMod
                         );
                     }

src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H

@@ -232,6 +232,7 @@ class addPatchCellLayer
             const label meshEdgeI,
             const label layerI,
             const label numEdgeFaces,
+            const labelList& meshFaces,
             polyTopoChange&
         ) const;
 

src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C

@@ -125,11 +125,11 @@ void Foam::combineFaces::regioniseFaces
 (
     const scalar minCos,
     const label cellI,
+    const labelList& cEdges,
     Map<label>& faceRegion
 ) const
 {
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-    const labelList& cEdges = mesh_.cellEdges()[cellI];
 
     forAll(cEdges, i)
     {
@@ -220,9 +220,10 @@ bool Foam::combineFaces::faceNeighboursValid
         return true;
     }
 
-    const labelListList& faceEdges = mesh_.faceEdges();
     const cell& cFaces = mesh_.cells()[cellI];
 
+    DynamicList<label> storage;
+
     // Test for face collapsing to edge since too many neighbours merged.
     forAll(cFaces, cFaceI)
     {
@@ -230,7 +231,7 @@ bool Foam::combineFaces::faceNeighboursValid
 
         if (!faceRegion.found(faceI))
         {
-            const labelList& fEdges = faceEdges[faceI];
+            const labelList& fEdges = mesh_.faceEdges(faceI, storage);
 
             // Count number of remaining faces neighbouring faceI. This has
             // to be 3 or more.
@@ -299,6 +300,8 @@ Foam::labelListList Foam::combineFaces::getMergeSets
 {
     // Lists of faces that can be merged.
     DynamicList<labelList> allFaceSets(boundaryCells.size() / 10);
+    // Storage for on-the-fly cell-edge addressing.
+    DynamicList<label> storage;
 
     // On all cells regionise the faces
     forAllConstIter(labelHashSet, boundaryCells, iter)
@@ -307,9 +310,11 @@ Foam::labelListList Foam::combineFaces::getMergeSets
 
         const cell& cFaces = mesh_.cells()[cellI];
 
+        const labelList& cEdges = mesh_.cellEdges(cellI, storage);
+
         // Region per face
         Map<label> faceRegion(cFaces.size());
-        regioniseFaces(featureCos, cellI, faceRegion);
+        regioniseFaces(featureCos, cellI, cEdges, faceRegion);
 
         // Now we have in faceRegion for every face the region with planar
         // face sharing the same region. We now check whether the resulting

src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.H

@@ -103,6 +103,7 @@ class combineFaces
         (
             const scalar minCos,
             const label cellI,
+            const labelList& cEdges,
             Map<label>& faceRegion
         ) const;
 

src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C

@@ -372,7 +372,7 @@ Foam::scalar Foam::hexRef8::getLevel0EdgeLength() const
         {
             const label cLevel = cellLevel_[cellI];
 
-            const labelList& cEdges = mesh_.cellEdges()[cellI];
+            const labelList& cEdges = mesh_.cellEdges(cellI);
 
             forAll(cEdges, i)
             {
@@ -447,7 +447,7 @@ Foam::scalar Foam::hexRef8::getLevel0EdgeLength() const
     {
         const label cLevel = cellLevel_[cellI];
 
-        const labelList& cEdges = mesh_.cellEdges()[cellI];
+        const labelList& cEdges = mesh_.cellEdges(cellI);
 
         forAll(cEdges, i)
         {
@@ -1190,6 +1190,10 @@ void Foam::hexRef8::createInternalFaces
     // From edge mid to face mids
     Map<edge> midPointToFaceMids(24);
 
+    // Storage for on-the-fly addressing
+    DynamicList<label> storage;
+
+
     // Running count of number of internal faces added so far.
     label nFacesAdded = 0;
 
@@ -1198,7 +1202,7 @@ void Foam::hexRef8::createInternalFaces
         label faceI = cFaces[i];
 
         const face& f = mesh_.faces()[faceI];
-        const labelList& fEdges = mesh_.faceEdges()[faceI];
+        const labelList& fEdges = mesh_.faceEdges(faceI, storage);
 
         // We are on the cellI side of face f. The face will have 1 or 4
         // cLevel points and lots of higher numbered ones.
@@ -1299,7 +1303,7 @@ void Foam::hexRef8::createInternalFaces
                     {
                         dumpCell(cellI);
 
-                        const labelList cPoints(cellPoints(cellI));
+                        const labelList& cPoints = mesh_.cellPoints(cellI);
 
                         FatalErrorIn("createInternalFaces(..)")
                             << "cell:" << cellI << " cLevel:" << cLevel
@@ -1372,7 +1376,7 @@ void Foam::hexRef8::createInternalFaces
                     {
                         dumpCell(cellI);
 
-                        const labelList cPoints(cellPoints(cellI));
+                        const labelList& cPoints = mesh_.cellPoints(cellI);
 
                         FatalErrorIn("createInternalFaces(..)")
                             << "cell:" << cellI << " cLevel:" << cLevel
@@ -1454,7 +1458,7 @@ void Foam::hexRef8::walkFaceToMid
 ) const
 {
     const face& f = mesh_.faces()[faceI];
-    const labelList& fEdges = mesh_.faceEdges()[faceI];
+    const labelList& fEdges = mesh_.faceEdges(faceI);
 
     label fp = startFp;
 
@@ -1503,7 +1507,7 @@ void Foam::hexRef8::walkFaceFromMid
 ) const
 {
     const face& f = mesh_.faces()[faceI];
-    const labelList& fEdges = mesh_.faceEdges()[faceI];
+    const labelList& fEdges = mesh_.faceEdges(faceI);
 
     label fp = f.rcIndex(startFp);
 
@@ -2013,27 +2017,6 @@ Foam::hexRef8::hexRef8
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-//- Get points of a cell (without using cellPoints addressing)
-Foam::labelList Foam::hexRef8::cellPoints(const label cellI) const
-{
-    // Pick up points of the cell
-    const cell& cFaces = mesh_.cells()[cellI];
-
-    labelHashSet cPoints(4*cFaces.size());
-
-    forAll(cFaces, i)
-    {
-        const face& f = mesh_.faces()[cFaces[i]];
-
-        forAll(f, fp)
-        {
-            cPoints.insert(f[fp]);
-        }
-    }
-    return cPoints.toc();
-}
-
-
 Foam::labelList Foam::hexRef8::consistentRefinement
 (
     const labelList& cellsToRefine,
@@ -2358,13 +2341,11 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement
         // as cell level purely for ease)
         labelList maxPointCount(mesh_.nPoints(), 0);
 
-        const labelListList& pointCells = mesh_.pointCells();
+        forAll(maxPointCount, pointI)
-
-        forAll(pointCells, pointI)
         {
             label& pLevel = maxPointCount[pointI];
 
-            const labelList& pCells = pointCells[pointI];
+            const labelList& pCells = mesh_.pointCells(pointI);
 
             forAll(pCells, i)
             {
@@ -2395,7 +2376,7 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement
             // Loop over all cells using the point and check whether their
             // refinement level is much less than the maximum.
 
-            const labelList& pCells = pointCells[pointI];
+            const labelList& pCells = mesh_.pointCells(pointI);
 
             forAll(pCells, pCellI)
             {
@@ -3121,7 +3102,7 @@ Foam::labelListList Foam::hexRef8::setRefinement
     {
         if (cellMidPoint[cellI] >= 0)
         {
-            const labelList& cEdges = mesh_.cellEdges()[cellI];
+            const labelList& cEdges = mesh_.cellEdges(cellI);
 
             forAll(cEdges, i)
             {
@@ -3458,7 +3439,7 @@ Foam::labelListList Foam::hexRef8::setRefinement
 
         forAll(pointLevel_, pointI)
         {
-            const labelList& pCells = mesh_.pointCells()[pointI];
+            const labelList& pCells = mesh_.pointCells(pointI);
 
             forAll(pCells, pCellI)
             {
@@ -3498,7 +3479,7 @@ Foam::labelListList Foam::hexRef8::setRefinement
                 {
                     dumpCell(cellI);
 
-                    const labelList cPoints(cellPoints(cellI));
+                    const labelList& cPoints = mesh_.cellPoints(cellI);
 
                     FatalErrorIn
                     (
@@ -3610,7 +3591,7 @@ Foam::labelListList Foam::hexRef8::setRefinement
         {
             if (edgeMidPoint[edgeI] >= 0)
             {
-                const labelList& eFaces = mesh_.edgeFaces()[edgeI];
+                const labelList& eFaces = mesh_.edgeFaces(edgeI);
 
                 forAll(eFaces, i)
                 {
@@ -3768,13 +3749,16 @@ Foam::labelListList Foam::hexRef8::setRefinement
             << endl;
     }
 
+    DynamicList<label> eFacesStorage;
+    DynamicList<label> fEdgesStorage;
+
     forAll(edgeMidPoint, edgeI)
     {
         if (edgeMidPoint[edgeI] >= 0)
         {
             // Split edge. Check that face not already handled above.
 
-            const labelList& eFaces = mesh_.edgeFaces()[edgeI];
+            const labelList& eFaces = mesh_.edgeFaces(edgeI, eFacesStorage);
 
             forAll(eFaces, i)
             {
@@ -3785,7 +3769,11 @@ Foam::labelListList Foam::hexRef8::setRefinement
                     // Unsplit face. Add edge splits to face.
 
                     const face& f = mesh_.faces()[faceI];
-                    const labelList& fEdges = mesh_.faceEdges()[faceI];
+                    const labelList& fEdges = mesh_.faceEdges
+                    (
+                        faceI,
+                        fEdgesStorage
+                    );
 
                     DynamicList<label> newFaceVerts(f.size());
 
@@ -4715,14 +4703,12 @@ void Foam::hexRef8::checkRefinementLevels
     // Check 2:1 across points (instead of faces)
     if (maxPointDiff != -1)
     {
-        const labelListList& pointCells = mesh_.pointCells();
-
         // Determine per point the max cell level.
         labelList maxPointLevel(mesh_.nPoints(), 0);
 
-        forAll(pointCells, pointI)
+        forAll(maxPointLevel, pointI)
         {
-            const labelList& pCells = pointCells[pointI];
+            const labelList& pCells = mesh_.pointCells(pointI);
 
             label& pLevel = maxPointLevel[pointI];
 
@@ -4747,7 +4733,7 @@ void Foam::hexRef8::checkRefinementLevels
         {
             label pointI = pointsToCheck[i];
 
-            const labelList& pCells = pointCells[pointI];
+            const labelList& pCells = mesh_.pointCells(pointI);
 
             forAll(pCells, i)
             {
@@ -4881,11 +4867,11 @@ Foam::labelList Foam::hexRef8::getSplitPoints() const
     labelList splitMasterLevel(mesh_.nPoints(), 0);
 
     // Unmark all with not 8 cells
-    const labelListList& pointCells = mesh_.pointCells();
+    //const labelListList& pointCells = mesh_.pointCells();
 
-    forAll(pointCells, pointI)
+    for (label pointI = 0; pointI < mesh_.nPoints(); pointI++)
     {
-        const labelList& pCells = pointCells[pointI];
+        const labelList& pCells = mesh_.pointCells(pointI);
 
         if (pCells.size() != 8)
         {
@@ -4898,8 +4884,7 @@ Foam::labelList Foam::hexRef8::getSplitPoints() const
 
     forAll(visibleCells, cellI)
     {
-        //const labelList& cPoints = mesh_.cellPoints()[cellI];
+        const labelList& cPoints = mesh_.cellPoints(cellI);
-        const labelList cPoints(cellPoints(cellI));
 
         if (visibleCells[cellI] != -1 && history_.parentIndex(cellI) >= 0)
         {
@@ -5104,7 +5089,7 @@ Foam::labelList Foam::hexRef8::consistentUnrefinement
         {
             if (unrefinePoint.get(pointI) == 1)
             {
-                const labelList& pCells = mesh_.pointCells()[pointI];
+                const labelList& pCells = mesh_.pointCells(pointI);
 
                 forAll(pCells, j)
                 {
@@ -5244,7 +5229,7 @@ Foam::labelList Foam::hexRef8::consistentUnrefinement
         {
             if (unrefinePoint.get(pointI) == 1)
             {
-                const labelList& pCells = mesh_.pointCells()[pointI];
+                const labelList& pCells = mesh_.pointCells(pointI);
 
                 forAll(pCells, j)
                 {
@@ -5329,7 +5314,7 @@ void Foam::hexRef8::setUnrefinement
 
         forAll(splitPointLabels, i)
         {
-            const labelList& pCells = mesh_.pointCells()[splitPointLabels[i]];
+            const labelList& pCells = mesh_.pointCells(splitPointLabels[i]);
 
             forAll(pCells, j)
             {
@@ -5395,7 +5380,7 @@ void Foam::hexRef8::setUnrefinement
 
         // Get original cell label
 
-        const labelList& pCells = mesh_.pointCells()[pointI];
+        const labelList& pCells = mesh_.pointCells(pointI);
 
         // Check
         if (pCells.size() != 8)
@@ -5463,7 +5448,7 @@ void Foam::hexRef8::setUnrefinement
     {
         label pointI = splitPointLabels[i];
 
-        const labelList& pCells = mesh_.pointCells()[pointI];
+        const labelList& pCells = mesh_.pointCells(pointI);
 
         label masterCellI = min(pCells);
 

src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.H

@@ -370,9 +370,6 @@ public:
 
         // Refinement
 
-            //- Helper:get points of a cell without using cellPoints addressing
-            labelList cellPoints(const label cellI) const;
-
             //- Given valid mesh and current cell level and proposed
             //  cells to refine calculate any clashes (due to 2:1) and return
             //  ok list of cells to refine.

src/dynamicMesh/polyTopoChange/polyTopoChange/removeFaces.C

@@ -83,6 +83,7 @@ Foam::label Foam::removeFaces::changeFaceRegion
     const labelList& nFacesPerEdge,
     const label faceI,
     const label newRegion,
+    const labelList& fEdges,
     labelList& faceRegion
 ) const
 {
@@ -94,27 +95,33 @@ Foam::label Foam::removeFaces::changeFaceRegion
 
         nChanged = 1;
 
+        // Storage for on-the-fly addressing
+        DynamicList<label> fe;
+        DynamicList<label> ef;
+
         // Step to neighbouring faces across edges that will get removed
-
-        const labelList& fEdges = mesh_.faceEdges()[faceI];
-
         forAll(fEdges, i)
         {
             label edgeI = fEdges[i];
 
             if (nFacesPerEdge[edgeI] >= 0 && nFacesPerEdge[edgeI] <= 2)
             {
-                const labelList& eFaces = mesh_.edgeFaces()[edgeI];
+                const labelList& eFaces = mesh_.edgeFaces(edgeI, ef);
 
                 forAll(eFaces, j)
                 {
+                    label nbrFaceI = eFaces[j];
+
+                    const labelList& fEdges1 = mesh_.faceEdges(nbrFaceI, fe);
+
                     nChanged += changeFaceRegion
                     (
                         cellRegion,
                         removedFace,
                         nFacesPerEdge,
-                        eFaces[j],
+                        nbrFaceI,
                         newRegion,
+                        fEdges1,
                         faceRegion
                     );
                 }
@@ -166,7 +173,7 @@ Foam::boolList Foam::removeFaces::getFacesAffected
     //  Mark faces affected by removal of edges
     forAllConstIter(labelHashSet, edgesToRemove, iter)
     {
-        const labelList& eFaces = mesh_.edgeFaces()[iter.key()];
+        const labelList& eFaces = mesh_.edgeFaces(iter.key());
 
         forAll(eFaces, eFaceI)
         {
@@ -814,6 +821,10 @@ void Foam::removeFaces::setRefinement
     // Number of connected face regions
     label nRegions = 0;
 
+    // Storage for on-the-fly addressing
+    DynamicList<label> fe;
+    DynamicList<label> ef;
+
 
     {
         const polyBoundaryMesh& patches = mesh_.boundaryMesh();
@@ -827,7 +838,7 @@ void Foam::removeFaces::setRefinement
         {
             label faceI = faceLabels[i];
 
-            const labelList& fEdges = mesh_.faceEdges()[faceI];
+            const labelList& fEdges = mesh_.faceEdges(faceI, fe);
 
             forAll(fEdges, i)
             {
@@ -835,8 +846,7 @@ void Foam::removeFaces::setRefinement
 
                 if (nFacesPerEdge[edgeI] == -1)
                 {
-                    nFacesPerEdge[edgeI] =
+                    nFacesPerEdge[edgeI] = mesh_.edgeFaces(edgeI, ef).size()-1;
-                        mesh_.edgeFaces()[edgeI].size()-1;
                 }
                 else
                 {
@@ -849,16 +859,15 @@ void Foam::removeFaces::setRefinement
         // Note that this only needs to be done for possibly coupled edges
         // so we could choose to loop only over boundary faces and use faceEdges
         // of those.
-        const labelListList& edgeFaces = mesh_.edgeFaces();
 
-        forAll(edgeFaces, edgeI)
+        forAll(mesh_.edges(), edgeI)
         {
             if (nFacesPerEdge[edgeI] == -1)
             {
                 // Edge not yet handled in loop above so is not used by any
                 // face to be removed.
 
-                const labelList& eFaces = edgeFaces[edgeI];
+                const labelList& eFaces = mesh_.edgeFaces(edgeI, ef);
 
                 if (eFaces.size() > 2)
                 {
@@ -922,7 +931,7 @@ void Foam::removeFaces::setRefinement
                 label f0 = -1;
                 label f1 = -1;
 
-                const labelList& eFaces = mesh_.edgeFaces()[edgeI];
+                const labelList& eFaces = mesh_.edgeFaces(edgeI, ef);
 
                 forAll(eFaces, i)
                 {
@@ -1152,6 +1161,7 @@ void Foam::removeFaces::setRefinement
                 nFacesPerEdge,
                 startFaceI,
                 nRegions,
+                mesh_.faceEdges(startFaceI, fe),
                 faceRegion
             );
 

src/dynamicMesh/polyTopoChange/polyTopoChange/removeFaces.H

@@ -93,6 +93,7 @@ class removeFaces
             const labelList& nFacesPerEdge,
             const label faceI,
             const label newRegion,
+            const labelList& fEdges,
             labelList& faceRegion
         ) const;
 

src/finiteVolume/Make/files

@@ -73,6 +73,7 @@ $(derivedFvPatchFields)/directMappedVelocityFluxFixedValue/directMappedVelocityF
 $(derivedFvPatchFields)/fan/fanFvPatchFields.C
 $(derivedFvPatchFields)/fixedFluxBuoyantPressure/fixedFluxBuoyantPressureFvPatchScalarField.C
 $(derivedFvPatchFields)/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.C
+$(derivedFvPatchFields)/fixedInternalValueFvPatchField/fixedInternalValueFvPatchFields.C
 $(derivedFvPatchFields)/fixedNormalSlip/fixedNormalSlipFvPatchFields.C
 $(derivedFvPatchFields)/fluxCorrectedVelocity/fluxCorrectedVelocityFvPatchVectorField.C
 $(derivedFvPatchFields)/freestream/freestreamFvPatchFields.C

src/finiteVolume/fields/fvPatchFields/derived/fixedInternalValueFvPatchField/fixedInternalValueFvPatchField.C

@@ -0,0 +1,105 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "fixedInternalValueFvPatchField.H"
+#include "fvPatchFieldMapper.H"
+#include "fvMatrix.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::fixedInternalValueFvPatchField<Type>::fixedInternalValueFvPatchField
+(
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF
+)
+:
+    zeroGradientFvPatchField<Type>(p, iF)
+{}
+
+
+template<class Type>
+Foam::fixedInternalValueFvPatchField<Type>::fixedInternalValueFvPatchField
+(
+    const fixedInternalValueFvPatchField<Type>& ptf,
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    zeroGradientFvPatchField<Type>(ptf, p, iF, mapper)
+{}
+
+
+template<class Type>
+Foam::fixedInternalValueFvPatchField<Type>::fixedInternalValueFvPatchField
+(
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    zeroGradientFvPatchField<Type>(p, iF, dict)
+{
+    fvPatchField<Type>::operator=(this->patchInternalField());
+}
+
+
+template<class Type>
+Foam::fixedInternalValueFvPatchField<Type>::fixedInternalValueFvPatchField
+(
+    const fixedInternalValueFvPatchField& fivpf
+)
+:
+    zeroGradientFvPatchField<Type>(fivpf)
+{}
+
+
+template<class Type>
+Foam::fixedInternalValueFvPatchField<Type>::fixedInternalValueFvPatchField
+(
+    const fixedInternalValueFvPatchField& fivpf,
+    const DimensionedField<Type, volMesh>& iF
+)
+:
+    zeroGradientFvPatchField<Type>(fivpf, iF)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::fixedInternalValueFvPatchField<Type>::manipulateMatrix
+(
+    fvMatrix<Type>& matrix
+)
+{
+    // Apply the patch internal field as a constraint in the matrix
+    matrix.setValues(this->patch().faceCells(), this->patchInternalField());
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/derived/fixedInternalValueFvPatchField/fixedInternalValueFvPatchField.H

@@ -0,0 +1,149 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::fixedInternalValueFvPatchField
+
+Description
+    Boundary condition providing mechanism to set boundary (cell) values
+    directly into a matrix, i.e. to set a constraint condition. Default
+    behaviour is to act as a zero gradient condition.
+
+SourceFiles
+    fixedInternalValueFvPatchField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fixedInternalValueFvPatchField_H
+#define fixedInternalValueFvPatchField_H
+
+#include "zeroGradientFvPatchField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+              Class fixedInternalValueFvPatchField Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class fixedInternalValueFvPatchField
+:
+    public zeroGradientFvPatchField<Type>
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("fixedInternalValue");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        fixedInternalValueFvPatchField
+        (
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        fixedInternalValueFvPatchField
+        (
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping the given fixedInternalValueFvPatchField<Type>
+        //  onto a new patch
+        fixedInternalValueFvPatchField
+        (
+            const fixedInternalValueFvPatchField<Type>&,
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        fixedInternalValueFvPatchField
+        (
+            const fixedInternalValueFvPatchField<Type>&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchField<Type> > clone() const
+        {
+            return tmp<fvPatchField<Type> >
+            (
+                new fixedInternalValueFvPatchField<Type>(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        fixedInternalValueFvPatchField
+        (
+            const fixedInternalValueFvPatchField<Type>&,
+            const DimensionedField<Type, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchField<Type> > clone
+        (
+            const DimensionedField<Type, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchField<Type> >
+            (
+                new fixedInternalValueFvPatchField<Type>(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Evaluation functions
+
+            //-Manipulate a matrix
+            virtual void manipulateMatrix(fvMatrix<Type>& matrix);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "fixedInternalValueFvPatchField.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/derived/fixedInternalValueFvPatchField/fixedInternalValueFvPatchFields.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,30 +22,24 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Global
-    wallOmega
-
-Description
-    Set wall dissipation in the omega matrix
-
 \*---------------------------------------------------------------------------*/
 
+#include "fixedInternalValueFvPatchFields.H"
+#include "fvPatchFields.H"
+#include "addToRunTimeSelectionTable.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
 {
-    const fvPatchList& patches = mesh_.boundary();
 
-    forAll(patches, patchi)
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-    {
-        const fvPatch& p = patches[patchi];
 
-        if (isType<wallFvPatch>(p))
+makePatchFields(fixedInternalValue);
-        {
+
-            omegaEqn().setValues
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-            (
+
-                p.faceCells(),
+} // End namespace Foam
-                omega_.boundaryField()[patchi].patchInternalField()
-            );
-        }
-    }
-}
 
 // ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/derived/fixedInternalValueFvPatchField/fixedInternalValueFvPatchFields.H

@@ -0,0 +1,50 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fixedInternalValueFvPatchFields_H
+#define fixedInternalValueFvPatchFields_H
+
+#include "fixedInternalValueFvPatchField.H"
+#include "fieldTypes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeFieldTypedefs(fixedInternalValue)
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.C

@@ -28,6 +28,7 @@ License
 #include "dictionary.H"
 #include "fvMesh.H"
 #include "fvPatchFieldMapper.H"
+//#include "fvMatrices.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -223,6 +224,13 @@ void Foam::fvPatchField<Type>::evaluate(const Pstream::commsTypes)
 }
 
 
+template<class Type>
+void Foam::fvPatchField<Type>::manipulateMatrix(fvMatrix<Type>& matrix)
+{
+    // do nothing
+}
+
+
 template<class Type>
 void Foam::fvPatchField<Type>::write(Ostream& os) const
 {

src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H

@@ -27,12 +27,12 @@ Class
 
 Description
     Abstract base class with a fat-interface to all derived classes
-    covering all possible ways in which they might be used. 
+    covering all possible ways in which they might be used.
 
     The first level of derivation is to basic patchFields which cover
     zero-gradient, fixed-gradient, fixed-value and mixed conditions.
 
-    The next level of derivation covers all the specialised typed with
+    The next level of derivation covers all the specialised types with
     specific evaluation proceedures, particularly with respect to specific
     fields.
 
@@ -66,6 +66,9 @@ class volMesh;
 template<class Type>
 class fvPatchField;
 
+template<class Type>
+class fvMatrix;
+
 template<class Type>
 Ostream& operator<<(Ostream&, const fvPatchField<Type>&);
 
@@ -407,6 +410,10 @@ public:
             }
 
 
+            //- Manipulate matrix
+            virtual void manipulateMatrix(fvMatrix<Type>& matrix);
+
+
         //- Write
         virtual void write(Ostream&) const;
 

src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C

@@ -592,6 +592,20 @@ void Foam::fvMatrix<Type>::relax()
 }
 
 
+template<class Type>
+void Foam::fvMatrix<Type>::boundaryManipulate
+(
+    typename GeometricField<Type, fvPatchField, volMesh>::
+        GeometricBoundaryField& bFields
+)
+{
+    forAll(bFields, patchI)
+    {
+        bFields[patchI].manipulateMatrix(*this);
+    }
+}
+
+
 template<class Type>
 Foam::tmp<Foam::scalarField> Foam::fvMatrix<Type>::D() const
 {

src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H

@@ -351,6 +351,13 @@ public:
             //  alpha is read from controlDict
             void relax();
 
+            //- Manipulate based on a boundary field
+            void boundaryManipulate
+            (
+                typename GeometricField<Type, fvPatchField, volMesh>::
+                    GeometricBoundaryField& values
+            );
+
             //- Construct and return the solver
             //  Solver controls read from Istream
             autoPtr<fvSolver> solver(Istream&);
@@ -387,7 +394,7 @@ public:
 
             //- Return the face-flux field from the matrix
             tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
-            flux() const;
+                flux() const;
 
 
     // Member operators

src/finiteVolume/fvMesh/wallDist/wallDistData.C

@@ -83,11 +83,11 @@ void Foam::wallDistData<TransferType>::correct()
     labelHashSet wallPatchIDs(getPatchIDs(wallPolyPatch::typeName));
 
     // Collect pointers to data on patches
-    List<Field<Type>*> patchData(mesh.boundaryMesh().size());
+    UPtrList<Field<Type> > patchData(mesh.boundaryMesh().size());
 
     forAll(field_.boundaryField(), patchI)
     {
-        patchData[patchI] = &(field_.boundaryField()[patchI]);
+        patchData.set(patchI, &field_.boundaryField()[patchI]);
     }
 
     // Do mesh wave

src/meshTools/cellDist/patchWave/patchDataWave.C

@@ -22,8 +22,6 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
 #include "patchDataWave.H"
@@ -50,8 +48,7 @@ void Foam::patchDataWave<TransferType>::setChangedFaces
         {
             const polyPatch& patch = mesh.boundaryMesh()[patchI];
 
-            const Field<Type>& patchField =
+            const Field<Type>& patchField = initialPatchValuePtrs_[patchI];
-                *initialPatchValuePtrs_[patchI];
 
             forAll(patch.faceCentres(), patchFaceI)
             {
@@ -176,7 +173,7 @@ Foam::patchDataWave<TransferType>::patchDataWave
 (
     const polyMesh& mesh,
     const labelHashSet& patchIDs,
-    const List<Field<Type>*>& initialPatchValuePtrs,
+    const UPtrList<Field<Type> >& initialPatchValuePtrs,
     const bool correctWalls
 )
 :

src/meshTools/cellDist/patchWave/patchDataWave.H

@@ -45,7 +45,7 @@ SourceFiles
 
 #include "cellDistFuncs.H"
 #include "FieldField.H"
-
+#include "UPtrList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -78,7 +78,7 @@ private:
         labelHashSet patchIDs_;
 
         //- Reference to initial extra data at patch faces
-        const List<Field<Type>*>& initialPatchValuePtrs_;
+        const UPtrList<Field<Type> >& initialPatchValuePtrs_;
 
         //- Do accurate distance calculation for near-wall cells.
         bool correctWalls_;
@@ -129,7 +129,7 @@ public:
         (
             const polyMesh& mesh,
             const labelHashSet& patchIDs,
-            const List<Field<Type>*>& initialPatchValuePtrs,
+            const UPtrList<Field<Type> >& initialPatchValuePtrs,
             bool correctWalls = true
         );
 

src/meshTools/sets/topoSets/faceSet.C

@@ -205,11 +205,12 @@ void faceSet::sync(const polyMesh& mesh)
 
     reduce(nAdded, sumOp<label>());
 
-    if (nAdded > 0)
+    //if (nAdded > 0)
-    {
+    //{
-        Info<< "Added an additional " << nAdded << " faces on coupled patches. "
+    //    Info<< "Added an additional " << nAdded
-            << "(processorPolyPatch, cyclicPolyPatch)" << endl;
+    //        << " faces on coupled patches. "
-    }
+    //        << "(processorPolyPatch, cyclicPolyPatch)" << endl;
+    //}
 }
 
 

src/postProcessing/Allwmake

@@ -6,5 +6,6 @@ wmake libo postCalc
 wmake libso forces
 wmake libso fieldAverage
 wmake libso foamCalcFunctions
+wmake libso minMaxFields
 
 # ----------------------------------------------------------------- end-of-file

src/postProcessing/minMaxFields/IOminMaxFields.H

@@ -22,30 +22,29 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Global
+Typedef
-    wallOmega
+    Foam::IOminMaxFields
 
 Description
-    Set wall dissipation in the omega matrix
+    Instance of the generic IOOutputFilter for minMaxFields.
 
 \*---------------------------------------------------------------------------*/
 
+#ifndef IOminMaxFields_H
+#define IOminMaxFields_H
+
+#include "minMaxFields.H"
+#include "IOOutputFilter.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
 {
-    const fvPatchList& patches = mesh_.boundary();
+    typedef IOOutputFilter<minMaxFields> IOminMaxFields;
-
-    forAll(patches, patchi)
-    {
-        const fvPatch& p = patches[patchi];
-
-        if (isType<wallFvPatch>(p))
-        {
-            omegaEqn().setValues
-            (
-                p.faceCells(),
-                omega_.boundaryField()[patchi].patchInternalField()
-            );
-        }
-    }
 }
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
 // ************************************************************************* //

src/postProcessing/minMaxFields/Make/files

@@ -0,0 +1,4 @@
+minMaxFields.C
+minMaxFieldsFunctionObject.C
+
+LIB = $(FOAM_LIBBIN)/libminMaxFields

src/postProcessing/minMaxFields/Make/options

@@ -0,0 +1,8 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+LIB_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools

src/postProcessing/minMaxFields/minMaxFields.C

@@ -0,0 +1,203 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "minMaxFields.H"
+#include "volFields.H"
+#include "dictionary.H"
+#include "Time.H"
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(minMaxFields, 0);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::minMaxFields::minMaxFields
+(
+    const word& name,
+    const objectRegistry& obr,
+    const dictionary& dict,
+    const bool loadFromFiles
+)
+:
+    name_(name),
+    obr_(obr),
+    active_(true),
+    log_(false),
+    fieldSet_(),
+    minMaxFieldsFilePtr_(NULL)
+{
+    // Check if the available mesh is an fvMesh otherise deactivate
+    if (!isA<fvMesh>(obr_))
+    {
+        active_ = false;
+        WarningIn
+        (
+            "minMaxFields::minMaxFields"
+            "(const objectRegistry& obr, const dictionary& dict)"
+        )   << "No fvMesh available, deactivating."
+            << endl;
+    }
+
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::minMaxFields::~minMaxFields()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::minMaxFields::read(const dictionary& dict)
+{
+    if (active_)
+    {
+        log_ = dict.lookupOrDefault<Switch>("log", false);
+
+        dict.lookup("fields") >> fieldSet_;
+    }
+}
+
+
+void Foam::minMaxFields::makeFile()
+{
+    // Create the minMaxFields file if not already created
+    if (!minMaxFieldsFilePtr_.valid())
+    {
+        if (debug)
+        {
+            Info<< "Creating minMaxFields file." << endl;
+        }
+
+        // File update
+        if (Pstream::master())
+        {
+            fileName minMaxFieldsDir;
+            if (Pstream::parRun())
+            {
+                // Put in undecomposed case (Note: gives problems for
+                // distributed data running)
+                minMaxFieldsDir =
+                    obr_.time().path()/".."/name_/obr_.time().timeName();
+            }
+            else
+            {
+                minMaxFieldsDir =
+                    obr_.time().path()/name_/obr_.time().timeName();
+            }
+
+            // Create directory if does not exist.
+            mkDir(minMaxFieldsDir);
+
+            // Open new file at start up
+            minMaxFieldsFilePtr_.reset
+            (
+                new OFstream(minMaxFieldsDir/(type() + ".dat"))
+            );
+
+            // Add headers to output data
+            writeFileHeader();
+        }
+    }
+}
+
+
+void Foam::minMaxFields::writeFileHeader()
+{
+    if (minMaxFieldsFilePtr_.valid())
+    {
+        minMaxFieldsFilePtr_()
+            << "# Time" << tab << "field" << tab << "min" << tab << "max"
+            << endl;
+    }
+}
+
+
+void Foam::minMaxFields::execute()
+{
+    // Do nothing - only valid on write
+}
+
+void Foam::minMaxFields::write()
+{
+    if (active_)
+    {
+        // Create the minMaxFields file if not already created
+        makeFile();
+
+        forAll(fieldSet_, fieldI)
+        {
+            calcMinMaxFields<scalar>(fieldSet_[fieldI]);
+            calcMinMaxFields<vector>(fieldSet_[fieldI]);
+            calcMinMaxFields<sphericalTensor>(fieldSet_[fieldI]);
+            calcMinMaxFields<symmTensor>(fieldSet_[fieldI]);
+            calcMinMaxFields<tensor>(fieldSet_[fieldI]);
+        }
+    }
+}
+
+
+template<>
+void Foam::minMaxFields::calcMinMaxFields<Foam::scalar>
+(
+    const word& fieldName
+)
+{
+    if (obr_.foundObject<volScalarField>(fieldName))
+    {
+        const scalarField& field = obr_.lookupObject<scalarField>(fieldName);
+        scalar minValue = min(field);
+        scalar maxValue = max(field);
+
+        reduce(minValue, minOp<scalar>());
+        reduce(maxValue, maxOp<scalar>());
+
+        if (Pstream::master())
+        {
+            minMaxFieldsFilePtr_() << obr_.time().value() << tab
+                << fieldName << tab << minValue << tab << maxValue << endl;
+
+            if (log_)
+            {
+                Info<< "minMaxFields output:" << nl
+                    << "    min(" << fieldName << ") = " << minValue << nl
+                    << "    max(" << fieldName << ") = " << maxValue << nl
+                    << endl;
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/postProcessing/minMaxFields/minMaxFields.H

@@ -0,0 +1,182 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::minMaxFields
+
+Description
+    Calculates scalar minimim and maximum field values.
+
+    For variables with rank > 0, computes the magnitude of the min/max
+    values.
+
+    Data written to the file \<timeDir\>/minMaxFields.dat
+
+
+SourceFiles
+    minMaxFields.C
+    IOminMaxFields.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef minMaxFields_H
+#define minMaxFields_H
+
+#include "primitiveFieldsFwd.H"
+#include "volFieldsFwd.H"
+#include "labelHashSet.H"
+#include "OFstream.H"
+#include "Switch.H"
+#include "pointFieldFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class objectRegistry;
+class dictionary;
+class mapPolyMesh;
+
+/*---------------------------------------------------------------------------*\
+                       Class minMaxFields Declaration
+\*---------------------------------------------------------------------------*/
+
+class minMaxFields
+{
+protected:
+
+    // Private data
+
+        //- Name of this set of forces,
+        //  Also used as the name of the probes directory.
+        word name_;
+
+        const objectRegistry& obr_;
+
+        //- on/off switch
+        bool active_;
+
+        //- Switch to send output to Info as well as to file
+        Switch log_;
+
+        //- Patches to integrate forces over
+        wordList fieldSet_;
+
+
+        //- Forces/moment file ptr
+        autoPtr<OFstream> minMaxFieldsFilePtr_;
+
+
+    // Private Member Functions
+
+        //- If the forces file has not been created create it
+        void makeFile();
+
+        //- Disallow default bitwise copy construct
+        minMaxFields(const minMaxFields&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const minMaxFields&);
+
+        //- Output file header information
+        virtual void writeFileHeader();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("minMaxFields");
+
+
+    // Constructors
+
+        //- Construct for given objectRegistry and dictionary.
+        //  Allow the possibility to load fields from files
+        minMaxFields
+        (
+            const word& name,
+            const objectRegistry&,
+            const dictionary&,
+            const bool loadFromFiles = false
+        );
+
+
+    // Destructor
+
+        virtual ~minMaxFields();
+
+
+    // Member Functions
+
+        //- Return name of the set of forces
+        virtual const word& name() const
+        {
+            return name_;
+        }
+
+        //- Read the forces data
+        virtual void read(const dictionary&);
+
+        //- Execute
+        virtual void execute();
+
+        //- Calculate the field min/max
+        template<class Type>
+        void calcMinMaxFields(const word& fieldName);
+
+        //- Write the minMaxFields
+        virtual void write();
+
+        //- Update for changes of mesh
+        virtual void updateMesh(const mapPolyMesh&)
+        {}
+
+        //- Update for changes of mesh
+        virtual void movePoints(const pointField&)
+        {}
+};
+
+
+// Template specialisation for scalar fields
+template<>
+void minMaxFields::calcMinMaxFields<scalar>(const word& fieldName);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "minMaxFieldsTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/postProcessing/minMaxFields/minMaxFieldsFunctionObject.C

@@ -24,15 +24,19 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "checkLock.H"
+#include "minMaxFieldsFunctionObject.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    globalCheckLock incompressibleTurbulenceModelsLock
+    defineNamedTemplateTypeNameAndDebug(minMaxFieldsFunctionObject, 0);
+
+    addToRunTimeSelectionTable
     (
-        "incompressibleTurbulenceModels"
+        functionObject,
+        minMaxFieldsFunctionObject,
+        dictionary
     );
 }
 

src/postProcessing/minMaxFields/minMaxFieldsFunctionObject.H

@@ -0,0 +1,55 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Typedef
+    Foam::minMaxFieldsFunctionObject
+
+Description
+    FunctionObject wrapper around minMaxFields to allow them to be created via
+    the functions list within controlDict.
+
+SourceFiles
+    minMaxFieldsFunctionObject.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef minMaxFieldsFunctionObject_H
+#define minMaxFieldsFunctionObject_H
+
+#include "minMaxFields.H"
+#include "OutputFilterFunctionObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef OutputFilterFunctionObject<minMaxFields>
+        minMaxFieldsFunctionObject;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/postProcessing/minMaxFields/minMaxFieldsTemplates.C

@@ -22,45 +22,40 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Global
-    kOmegaWallViscosity
-
-Description
-    Calculate wall viscosity from wall-functions.
-
 \*---------------------------------------------------------------------------*/
 
+#include "minMaxFields.H"
+#include "volFields.H"
+#include "dictionary.H"
+#include "Time.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::minMaxFields::calcMinMaxFields(const word& fieldName)
 {
-    scalar Cmu25 = pow(Cmu_.value(), 0.25);
+    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
 
-    const fvPatchList& patches = mesh_.boundary();
+    if (obr_.foundObject<fieldType>(fieldName))
-
-    forAll(patches, patchi)
     {
-        const fvPatch& curPatch = patches[patchi];
+        const fieldType& field = obr_.lookupObject<fieldType>(fieldName);
+        scalar minValue = min(mag(field)).value();
+        scalar maxValue = max(mag(field)).value();
 
-        if (isType<wallFvPatch>(curPatch))
+        reduce(minValue, minOp<scalar>());
+        reduce(maxValue, maxOp<scalar>());
+
+        if (Pstream::master())
         {
-            const scalarField& nuw = nu().boundaryField()[patchi];
+            minMaxFieldsFilePtr_() << obr_.time().value() << tab
-            scalarField& nutw = nut_.boundaryField()[patchi];
+                << fieldName << tab << minValue << tab << maxValue << endl;
 
-            forAll(curPatch, facei)
+            if (log_)
             {
-                label faceCelli = curPatch.faceCells()[facei];
+                Info<< "minMaxFields output:" << nl
-
+                    << "    min(mag(" << fieldName << ")) = " << minValue << nl
-                scalar yPlus =
+                    << "    max(mag(" << fieldName << ")) = " << maxValue << nl
-                    Cmu25*y_[patchi][facei]*sqrt(k_[faceCelli])/nuw[facei];
+                    << endl;
-
-                if (yPlus > yPlusLam_)
-                {
-                    nutw[facei] =
-                         nuw[facei]
-                        *(yPlus*kappa_.value()/log(E_.value()*yPlus) - 1);
-                }
-                else
-                {
-                    nutw[facei] = 0.0;
-                }
             }
         }
     }

src/turbulenceModels/LES/compressible/Make/files

@@ -1,5 +1,3 @@
-wallFunc/muSgsWallFunction/muSgsWallFunctionFvPatchScalarField.C
-
 LESModel/LESModel.C
 LESModel/newLESModel.C
 GenEddyVisc/GenEddyVisc.C
@@ -12,4 +10,10 @@ dynOneEqEddy/dynOneEqEddy.C
 DeardorffDiffStress/DeardorffDiffStress.C
 SpalartAllmaras/SpalartAllmaras.C
 
+/* Wall functions */
+wallFunctions=derivedFvPatchFields/wallFunctions
+
+muSgsWallFunctions=$(wallFunctions)/muSgsWallFunctions
+$(muSgsWallFunctions)/muSgsWallFunction/muSgsWallFunctionFvPatchScalarField.C
+
 LIB = $(FOAM_LIBBIN)/libcompressibleLESModels

src/turbulenceModels/LES/incompressible/Make/files

@@ -1,5 +1,4 @@
 vanDriestDelta/vanDriestDelta.C
-wallFunc/nuSgsWallFunction/nuSgsWallFunctionFvPatchScalarField.C
 
 LESModel/LESModel.C
 LESModel/newLESModel.C
@@ -25,4 +24,12 @@ dynMixedSmagorinsky/dynMixedSmagorinsky.C
 
 /*Smagorinsky2/Smagorinsky2.C*/
 
+
+/* Wall functions */
+wallFunctions=derivedFvPatchFields/wallFunctions
+
+nuSgsWallFunctions=$(wallFunctions)/nuSgsWallFunctions
+$(nuSgsWallFunctions)/nuSgsWallFunction/nuSgsWallFunctionFvPatchScalarField.C
+
+
 LIB = $(FOAM_LIBBIN)/libincompressibleLESModels

src/turbulenceModels/RAS/compressible/LRR/LRR.C

@@ -28,6 +28,8 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "wallFvPatch.H"
 
+#include "backwardsCompatibilityWallFunctions.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -165,9 +167,8 @@ LRR::LRR
             IOobject::MUST_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh_
+        autoCreateR("R", mesh_)
     ),
-
     k_
     (
         IOobject
@@ -175,12 +176,11 @@ LRR::LRR
             "k",
             runTime_.timeName(),
             mesh_,
-            IOobject::MUST_READ,
+            IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh_
+        autoCreateK("k", mesh_)
     ),
-
     epsilon_
     (
         IOobject
@@ -188,12 +188,11 @@ LRR::LRR
             "epsilon",
             runTime_.timeName(),
             mesh_,
-            IOobject::MUST_READ,
+            IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh_
+        autoCreateEpsilon("epsilon", mesh_)
     ),
-
     mut_
     (
         IOobject
@@ -202,12 +201,13 @@ LRR::LRR
             runTime_.timeName(),
             mesh_,
             IOobject::NO_READ,
-            IOobject::NO_WRITE
+            IOobject::AUTO_WRITE
         ),
-        Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_)
+        autoCreateMut("mut", mesh_)
     )
 {
-#   include "wallViscosityI.H"
+    mut_ = Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    mut_.correctBoundaryConditions();
 
     if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
     {
@@ -311,15 +311,17 @@ void LRR::correct()
     {
         // Re-calculate viscosity
         mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
+        mut_.correctBoundaryConditions();
         return;
     }
 
     RASModel::correct();
 
     volSymmTensorField P = -twoSymm(R_ & fvc::grad(U_));
-    volScalarField G = 0.5*tr(P);
+    volScalarField G("G", 0.5*tr(P));
 
-#   include "wallFunctionsI.H"
+    // Update espsilon and G at the wall
+    epsilon_.boundaryField().updateCoeffs();
 
     // Dissipation equation
     tmp<fvScalarMatrix> epsEqn
@@ -335,7 +337,7 @@ void LRR::correct()
 
     epsEqn().relax();
 
-#   include "wallDissipationI.H"
+    epsEqn().boundaryManipulate(epsilon_.boundaryField());
 
     solve(epsEqn);
     bound(epsilon_, epsilon0_);
@@ -398,8 +400,7 @@ void LRR::correct()
 
     // Re-calculate viscosity
     mut_ = rho_*Cmu_*sqr(k_)/epsilon_;
-
+    mut_.correctBoundaryConditions();
-#   include "wallViscosityI.H"
 
 
     // Correct wall shear stresses

src/turbulenceModels/RAS/compressible/LaunderGibsonRSTM/LaunderGibsonRSTM.C

@@ -30,6 +30,8 @@ License
 #include "wallDist.H"
 #include "wallDistReflection.H"
 
+#include "backwardsCompatibilityWallFunctions.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -187,9 +189,8 @@ LaunderGibsonRSTM::LaunderGibsonRSTM
             IOobject::MUST_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh_
+        autoCreateR("R", mesh_)
     ),
-
     k_
     (
         IOobject
@@ -197,12 +198,11 @@ LaunderGibsonRSTM::LaunderGibsonRSTM
             "k",
             runTime_.timeName(),
             mesh_,
-            IOobject::MUST_READ,
+            IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh_
+        autoCreateK("k", mesh_)
     ),
-
     epsilon_
     (
         IOobject
@@ -210,12 +210,11 @@ LaunderGibsonRSTM::LaunderGibsonRSTM
             "epsilon",
             runTime_.timeName(),
             mesh_,
-            IOobject::MUST_READ,
+            IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh_
+        autoCreateEpsilon("epsilon", mesh_)
     ),
-
     mut_
     (
         IOobject
@@ -224,12 +223,13 @@ LaunderGibsonRSTM::LaunderGibsonRSTM
             runTime_.timeName(),
             mesh_,
             IOobject::NO_READ,
-            IOobject::NO_WRITE
+            IOobject::AUTO_WRITE
         ),
-        Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_)
+        autoCreateMut("mut", mesh_)
     )
 {
-#   include "wallViscosityI.H"
+    mut_ = Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    mut_.correctBoundaryConditions();
 
     if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
     {
@@ -336,6 +336,7 @@ void LaunderGibsonRSTM::correct()
     {
         // Re-calculate viscosity
         mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
+        mut_.correctBoundaryConditions();
         return;
     }
 
@@ -347,9 +348,10 @@ void LaunderGibsonRSTM::correct()
     }
 
     volSymmTensorField P = -twoSymm(R_ & fvc::grad(U_));
-    volScalarField G = 0.5*tr(P);
+    volScalarField G("G", 0.5*tr(P));
 
-#   include "wallFunctionsI.H"
+    // Update espsilon and G at the wall
+    epsilon_.boundaryField().updateCoeffs();
 
     // Dissipation equation
     tmp<fvScalarMatrix> epsEqn
@@ -365,7 +367,7 @@ void LaunderGibsonRSTM::correct()
 
     epsEqn().relax();
 
-#   include "wallDissipationI.H"
+    epsEqn().boundaryManipulate(epsilon_.boundaryField());
 
     solve(epsEqn);
     bound(epsilon_, epsilon0_);
@@ -437,9 +439,7 @@ void LaunderGibsonRSTM::correct()
 
     // Re-calculate turbulent viscosity
     mut_ = Cmu_*rho_*sqr(k_)/epsilon_;
-
+    mut_.correctBoundaryConditions();
-
-#   include "wallViscosityI.H"
 
 
     // Correct wall shear stresses

src/turbulenceModels/RAS/compressible/Make/files

@@ -12,11 +12,29 @@ SpalartAllmaras/SpalartAllmaras.C
 kOmegaSST/kOmegaSST.C
 
 /* Wall functions */
-wallFunctions/mutWallFunctions/mutStandardRoughWallFunction/mutStandardRoughWallFunctionFvPatchScalarField.C
+wallFunctions = derivedFvPatchFields/wallFunctions
+
+mutWallFunctions = $(wallFunctions)/mutWallFunctions
+$(mutWallFunctions)/mutWallFunction/mutWallFunctionFvPatchScalarField.C
+$(mutWallFunctions)/mutRoughWallFunction/mutRoughWallFunctionFvPatchScalarField.C
+$(mutWallFunctions)/mutSpalartAllmarasWallFunction/mutSpalartAllmarasWallFunctionFvPatchScalarField.C
+$(mutWallFunctions)/mutSpalartAllmarasStandardWallFunction/mutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C
+$(mutWallFunctions)/mutSpalartAllmarasStandardRoughWallFunction/mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C
+
+epsilonWallFunctions = $(wallFunctions)/epsilonWallFunctions
+$(epsilonWallFunctions)/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C
+
+omegaWallFunctions = $(wallFunctions)/omegaWallFunctions
+$(omegaWallFunctions)/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
+
+kQRWallFunctions = $(wallFunctions)/kQRWallFunctions
+$(kQRWallFunctions)/kQRWallFunction/kQRWallFunctionFvPatchFields.C
 
 /* Patch fields */
 derivedFvPatchFields/turbulentHeatFluxTemperature/turbulentHeatFluxTemperatureFvPatchScalarField.C
 derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.C
 derivedFvPatchFields/turbulentMixingLengthFrequencyInlet/turbulentMixingLengthFrequencyInletFvPatchScalarField.C
 
+backwardsCompatibilityWallFunctions/backwardsCompatibilityWallFunctions.C
+
 LIB = $(FOAM_LIBBIN)/libcompressibleRASModels

src/turbulenceModels/RAS/compressible/RASModel/RASModel.C

@@ -87,6 +87,7 @@ RASModel::RASModel
     printCoeffs_(lookupOrDefault<Switch>("printCoeffs", false)),
     coeffDict_(subDict(type + "Coeffs")),
 
+    wallFunctionDict_(subDict("wallFunctionCoeffs")),
     kappa_
     (
         dimensioned<scalar>::lookupOrAddToDict
@@ -105,6 +106,15 @@ RASModel::RASModel
             9.0
         )
     ),
+    Cmu_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cmu",
+            wallFunctionDict_,
+            0.09
+        )
+    ),
 
     yPlusLam_(yPlusLam(kappa_.value(), E_.value())),
 
@@ -118,7 +128,7 @@ RASModel::RASModel
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-scalar RASModel::yPlusLam(const scalar kappa, const scalar E)
+scalar RASModel::yPlusLam(const scalar kappa, const scalar E) const
 {
     scalar ypl = 11.0;
 

src/turbulenceModels/RAS/compressible/RASModel/RASModel.H

@@ -91,10 +91,11 @@ protected:
         Switch printCoeffs_;
         dictionary coeffDict_;
 
+        dictionary wallFunctionDict_;
         dimensionedScalar kappa_;
         dimensionedScalar E_;
+        dimensionedScalar Cmu_;
 
-        scalar yPlusLam(const scalar kappa, const scalar E);
         scalar yPlusLam_;
 
         dimensionedScalar k0_;
@@ -237,6 +238,21 @@ public:
                 return E_;
             }
 
+            //- Return Cmu for use in wall-functions
+            dimensionedScalar Cmu() const
+            {
+                return Cmu_;
+            }
+
+            //- Return the near wall distances
+            const nearWallDist& y() const
+            {
+                return y_;
+            }
+
+            //- Calculate y+ at the edge of the laminar sublayer
+            scalar yPlusLam(const scalar kappa, const scalar E) const;
+
             //- Return y+ at the edge of the laminar sublayer
             //  for use in wall-functions
             scalar yPlusLam() const
@@ -250,6 +266,12 @@ public:
                 return coeffDict_;
             }
 
+            //- Const access to the wall functions coefficients dictionary
+            const dictionary& walLFunctionDict() const
+            {
+                return wallFunctionDict_;
+            }
+
 
         //- Return the laminar viscosity
         const volScalarField& mu() const

src/turbulenceModels/RAS/compressible/RNGkEpsilon/RNGkEpsilon.C

@@ -28,6 +28,8 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "wallFvPatch.H"
 
+#include "backwardsCompatibilityWallFunctions.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -144,12 +146,11 @@ RNGkEpsilon::RNGkEpsilon
             "k",
             runTime_.timeName(),
             mesh_,
-            IOobject::MUST_READ,
+            IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh_
+        autoCreateK("k", mesh_)
     ),
-
     epsilon_
     (
         IOobject
@@ -157,12 +158,11 @@ RNGkEpsilon::RNGkEpsilon
             "epsilon",
             runTime_.timeName(),
             mesh_,
-            IOobject::MUST_READ,
+            IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh_
+        autoCreateEpsilon("epsilon", mesh_)
     ),
-
     mut_
     (
         IOobject
@@ -171,12 +171,13 @@ RNGkEpsilon::RNGkEpsilon
             runTime_.timeName(),
             mesh_,
             IOobject::NO_READ,
-            IOobject::NO_WRITE
+            IOobject::AUTO_WRITE
         ),
-        Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_)
+        autoCreateMut("mut", mesh_)
     )
 {
-#   include "wallViscosityI.H"
+    mut_ = Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    mut_.correctBoundaryConditions();
 
     printCoeffs();
 }
@@ -263,6 +264,7 @@ void RNGkEpsilon::correct()
     {
         // Re-calculate viscosity
         mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
+        mut_.correctBoundaryConditions();
         return;
     }
 
@@ -279,7 +281,7 @@ void RNGkEpsilon::correct()
     volScalarField S2 = (tgradU() && dev(twoSymm(tgradU())));
     tgradU.clear();
 
-    volScalarField G = mut_*S2;
+    volScalarField G("G", mut_*S2);
 
     volScalarField eta = sqrt(mag(S2))*k_/epsilon_;
     volScalarField eta3 = eta*sqr(eta);
@@ -287,7 +289,8 @@ void RNGkEpsilon::correct()
     volScalarField R =
         ((eta*(-eta/eta0_ + scalar(1)))/(beta_*eta3 + scalar(1)));
 
-#   include "wallFunctionsI.H"
+    // Update espsilon and G at the wall
+    epsilon_.boundaryField().updateCoeffs();
 
     // Dissipation equation
     tmp<fvScalarMatrix> epsEqn
@@ -303,7 +306,7 @@ void RNGkEpsilon::correct()
 
     epsEqn().relax();
 
-#   include "wallDissipationI.H"
+    epsEqn().boundaryManipulate(epsilon_.boundaryField());
 
     solve(epsEqn);
     bound(epsilon_, epsilon0_);
@@ -328,9 +331,7 @@ void RNGkEpsilon::correct()
 
     // Re-calculate viscosity
     mut_ = rho_*Cmu_*sqr(k_)/epsilon_;
-
+    mut_.correctBoundaryConditions();
-#   include "wallViscosityI.H"
-
 }
 
 

src/turbulenceModels/RAS/compressible/backwardsCompatibilityWallFunctions/backwardsCompatibilityWallFunctions.C

@@ -0,0 +1,215 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "backwardsCompatibilityWallFunctions.H"
+
+#include "calculatedFvPatchField.H"
+#include "mutWallFunctionFvPatchScalarField.H"
+#include "epsilonWallFunctionFvPatchScalarField.H"
+#include "kQRWallFunctionFvPatchField.H"
+#include "omegaWallFunctionFvPatchScalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+tmp<volScalarField> autoCreateMut
+(
+    const word& fieldName,
+    const fvMesh& mesh
+)
+{
+    IOobject mutHeader
+    (
+        fieldName,
+        mesh.time().timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE,
+        false
+    );
+
+    if (mutHeader.headerOk())
+    {
+        return tmp<volScalarField>(new volScalarField(mutHeader, mesh));
+    }
+    else
+    {
+        Info<< "--> Upgrading " << fieldName << " to employ run-time "
+            << "selectable wall functions" << endl;
+
+        const fvBoundaryMesh& bm = mesh.boundary();
+
+        wordList mutBoundaryTypes(bm.size());
+
+        forAll(bm, patchI)
+        {
+            if (isType<wallFvPatch>(bm[patchI]))
+            {
+                mutBoundaryTypes[patchI] =
+                    RASModels::mutWallFunctionFvPatchScalarField::typeName;
+            }
+            else
+            {
+                mutBoundaryTypes[patchI] =
+                    calculatedFvPatchField<scalar>::typeName;
+            }
+        }
+
+        tmp<volScalarField> mut
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    fieldName,
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                mesh,
+                dimensionedScalar("zero", dimDensity*dimArea/dimTime, 0.0),
+                mutBoundaryTypes
+            )
+        );
+
+        Info<< "    Writing updated " << fieldName << endl;
+        mut().write();
+
+        return mut;
+    }
+}
+
+
+tmp<volScalarField> autoCreateEpsilon
+(
+    const word& fieldName,
+    const fvMesh& mesh
+)
+{
+    return
+        autoCreateWallFunctionField
+        <
+            scalar,
+            RASModels::epsilonWallFunctionFvPatchScalarField
+        >
+        (
+            fieldName,
+            mesh
+        );
+}
+
+
+tmp<volScalarField> autoCreateOmega
+(
+    const word& fieldName,
+    const fvMesh& mesh
+)
+{
+    return
+        autoCreateWallFunctionField
+        <
+            scalar,
+            RASModels::omegaWallFunctionFvPatchScalarField
+        >
+        (
+            fieldName,
+            mesh
+        );
+}
+
+
+tmp<volScalarField> autoCreateK
+(
+    const word& fieldName,
+    const fvMesh& mesh
+)
+{
+    return
+        autoCreateWallFunctionField
+        <
+            scalar,
+            RASModels::kQRWallFunctionFvPatchField<scalar>
+        >
+        (
+            fieldName,
+            mesh
+        );
+}
+
+
+tmp<volScalarField> autoCreateQ
+(
+    const word& fieldName,
+    const fvMesh& mesh
+)
+{
+    return
+        autoCreateWallFunctionField
+        <
+            scalar,
+            RASModels::kQRWallFunctionFvPatchField<scalar>
+        >
+        (
+            fieldName,
+            mesh
+        );
+}
+
+
+tmp<volSymmTensorField> autoCreateR
+(
+    const word& fieldName,
+    const fvMesh& mesh
+)
+{
+    return
+        autoCreateWallFunctionField
+        <
+            symmTensor,
+            RASModels::kQRWallFunctionFvPatchField<symmTensor>
+        >
+        (
+            fieldName,
+            mesh
+        );
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace compressible
+} // End namespace Foam
+
+// ************************************************************************* //
+

src/turbulenceModels/RAS/compressible/backwardsCompatibilityWallFunctions/backwardsCompatibilityWallFunctions.H

@@ -0,0 +1,116 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::compressible
+
+Description
+    Auto creation of fields to provide backwards compatibility with
+    runtime selectable wall functions
+
+SourceFiles
+    backwardsCompatibilityWallFunctions.C
+    backwardsCompatibilityWallFunctionsTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef backwardsCompatibilityWallFunctions_H
+#define backwardsCompatibilityWallFunctions_H
+
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+    //- mut
+    tmp<volScalarField> autoCreateMut
+    (
+        const word& fieldName,
+        const fvMesh& mesh
+    );
+
+    //- epsilon
+    tmp<volScalarField> autoCreateEpsilon
+    (
+        const word& fieldName,
+        const fvMesh& mesh
+    );
+
+    //- omega
+    tmp<volScalarField> autoCreateOmega
+    (
+        const word& fieldName,
+        const fvMesh& mesh
+    );
+
+    //- k
+    tmp<volScalarField> autoCreateK
+    (
+        const word& fieldName,
+        const fvMesh& mesh
+    );
+
+    //- Q
+    tmp<volScalarField> autoCreateQ
+    (
+        const word& fieldName,
+        const fvMesh& mesh
+    );
+
+    //- R
+    tmp<volSymmTensorField> autoCreateR
+    (
+        const word& fieldName,
+        const fvMesh& mesh
+    );
+
+    //- Helper function to create the new field
+    template<class Type, class PatchType>
+    tmp<GeometricField<Type, fvPatchField, volMesh> >
+    autoCreateWallFunctionField
+    (
+        const word& fieldName,
+        const fvMesh& mesh
+    );
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace compressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "backwardsCompatibilityWallFunctionsTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/turbulenceModels/RAS/compressible/backwardsCompatibilityWallFunctions/backwardsCompatibilityWallFunctionsTemplates.C

@@ -0,0 +1,166 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "backwardsCompatibilityWallFunctions.H"
+#include "Time.H"
+
+#include "wallPolyPatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type, class PatchType>
+tmp<GeometricField<Type, fvPatchField, volMesh> >
+autoCreateWallFunctionField
+(
+    const word& fieldName,
+    const fvMesh& mesh
+)
+{
+    IOobject mutHeader
+    (
+        "mut",
+        mesh.time().timeName(),
+        mesh,
+        IOobject::MUST_READ
+    );
+
+    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+
+    if (mutHeader.headerOk())
+    {
+        return tmp<fieldType>
+        (
+            new fieldType
+            (
+                IOobject
+                (
+                    fieldName,
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                mesh
+            )
+        );
+    }
+    else
+    {
+        Info<< "--> Upgrading " << fieldName << " to employ run-time "
+            << "selectable wall functions" << endl;
+
+        // Read existing epsilon field
+        tmp<fieldType> fieldOrig
+        (
+            new fieldType
+            (
+                IOobject
+                (
+                    fieldName,
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                mesh
+            )
+        );
+
+        PtrList<fvPatchField<Type> > newPatchFields(mesh.boundary().size());
+
+        forAll(newPatchFields, patchI)
+        {
+            if (isType<wallPolyPatch>(mesh.boundaryMesh()[patchI]))
+            {
+                newPatchFields.set
+                (
+                    patchI,
+                    new PatchType
+                    (
+                        mesh.boundary()[patchI],
+                        fieldOrig().dimensionedInternalField()
+                    )
+                );
+                newPatchFields[patchI] == fieldOrig().boundaryField()[patchI];
+            }
+            else
+            {
+                newPatchFields.set
+                (
+                    patchI,
+                    fieldOrig().boundaryField()[patchI].clone()
+                );
+            }
+        }
+
+        tmp<fieldType> fieldNew
+        (
+            new fieldType
+            (
+                IOobject
+                (
+                    fieldName,
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                mesh,
+                fieldOrig().dimensions(),
+                fieldOrig().internalField(),
+                newPatchFields
+            )
+        );
+
+        Info<< "    Writing backup of original " << fieldName << " to "
+            << fieldName << ".old" << endl;
+        fieldOrig().rename(fieldName + ".old");
+        fieldOrig().write();
+
+        Info<< "    Writing updated " << fieldName << endl;
+        fieldNew().write();
+
+        return fieldNew;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace compressible
+} // End namespace Foam
+
+// ************************************************************************* //

src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C

@@ -0,0 +1,218 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "epsilonWallFunctionFvPatchScalarField.H"
+#include "RASModel.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+#include "wallFvPatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void epsilonWallFunctionFvPatchScalarField::checkType()
+{
+    if (!isA<wallFvPatch>(patch()))
+    {
+        FatalErrorIn("epsilonWallFunctionFvPatchScalarField::checkType()")
+            << "Invalid wall function specification" << nl
+            << "    Patch type for patch " << patch().name()
+            << " must be wall" << nl
+            << "    Current patch type is " << patch().type() << nl << endl
+            << abort(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedInternalValueFvPatchField<scalar>(p, iF)
+{
+    checkType();
+}
+
+
+epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
+(
+    const epsilonWallFunctionFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedInternalValueFvPatchField<scalar>(ptf, p, iF, mapper)
+{
+    checkType();
+}
+
+
+epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedInternalValueFvPatchField<scalar>(p, iF, dict)
+{
+    checkType();
+}
+
+
+epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
+(
+    const epsilonWallFunctionFvPatchScalarField& ewfpsf
+)
+:
+    fixedInternalValueFvPatchField<scalar>(ewfpsf)
+{
+    checkType();
+}
+
+
+epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
+(
+    const epsilonWallFunctionFvPatchScalarField& ewfpsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedInternalValueFvPatchField<scalar>(ewfpsf, iF)
+{
+    checkType();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void epsilonWallFunctionFvPatchScalarField::updateCoeffs()
+{
+    const RASModel& ras = db().lookupObject<RASModel>("RASProperties");
+
+    const scalar Cmu = ras.Cmu().value();
+    const scalar Cmu25 = pow(Cmu, 0.25);
+    const scalar Cmu75 = pow(Cmu, 0.75);
+    const scalar kappa = ras.kappa().value();
+    const scalar yPlusLam = ras.yPlusLam();
+
+    const scalarField& y = ras.y()[patch().index()];
+
+    volScalarField& G = const_cast<volScalarField&>
+        (db().lookupObject<volScalarField>("G"));
+
+    volScalarField& epsilon = const_cast<volScalarField&>
+        (db().lookupObject<volScalarField>("epsilon"));
+
+    const volScalarField& k = db().lookupObject<volScalarField>("k");
+
+    const scalarField& rhow =
+        patch().lookupPatchField<volScalarField, scalar>("rho");
+
+    const scalarField& muw =
+        patch().lookupPatchField<volScalarField, scalar>("mu");
+
+    const scalarField& mutw =
+        patch().lookupPatchField<volScalarField, scalar>("mut");
+
+    const fvPatchVectorField& Uw =
+        patch().lookupPatchField<volVectorField, vector>("U");
+
+    const scalarField magGradUw = mag(Uw.snGrad());
+
+    // Set epsilon and G
+    forAll(mutw, faceI)
+    {
+        label faceCellI = patch().faceCells()[faceI];
+
+        scalar yPlus =
+            Cmu25*y[faceI]*sqrt(k[faceCellI])
+           /(muw[faceI]/rhow[faceI]);
+
+        epsilon[faceCellI] = Cmu75*pow(k[faceCellI], 1.5)/(kappa*y[faceI]);
+
+        if (yPlus > yPlusLam)
+        {
+            G[faceCellI] =
+                (mutw[faceI] + muw[faceI])
+               *magGradUw[faceI]
+               *Cmu25*sqrt(k[faceCellI])
+               /(kappa*y[faceI]);
+        }
+        else
+        {
+            G[faceCellI] = 0.0;
+        }
+    }
+
+    // TODO: perform averaging for cells sharing more than one boundary face
+}
+
+
+void epsilonWallFunctionFvPatchScalarField::evaluate
+(
+    const Pstream::commsTypes commsType
+)
+{
+    fixedInternalValueFvPatchField<scalar>::evaluate(commsType);
+}
+
+
+void epsilonWallFunctionFvPatchScalarField::write(Ostream& os) const
+{
+    fixedInternalValueFvPatchField<scalar>::write(os);
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField
+(
+    fvPatchScalarField,
+    epsilonWallFunctionFvPatchScalarField
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// ************************************************************************* //

src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.H

@@ -0,0 +1,164 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::compressible::RASModels::epsilonWallFunctionFvPatchScalarField
+
+Description
+    Boundary condition for epsilon when using wall functions
+    - calculates epsilon and G
+    - epsilon values added directly into the matrix to act as a constraint
+
+SourceFiles
+    epsilonWallFunctionFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef epsilonWallFunctionFvPatchScalarField_H
+#define epsilonWallFunctionFvPatchScalarField_H
+
+#include "fixedInternalValueFvPatchField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+           Class epsilonWallFunctionFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class epsilonWallFunctionFvPatchScalarField
+:
+    public fixedInternalValueFvPatchField<scalar>
+{
+
+    // Private member functions
+
+        //- Check the type of the patch
+        void checkType();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("epsilonWallFunction");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        epsilonWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        epsilonWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  epsilonWallFunctionFvPatchScalarField
+        //  onto a new patch
+        epsilonWallFunctionFvPatchScalarField
+        (
+            const epsilonWallFunctionFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        epsilonWallFunctionFvPatchScalarField
+        (
+            const epsilonWallFunctionFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new epsilonWallFunctionFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        epsilonWallFunctionFvPatchScalarField
+        (
+            const epsilonWallFunctionFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new epsilonWallFunctionFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+            //- Evaluate the patchField
+            virtual void evaluate(const Pstream::commsTypes);
+
+
+        // I-O
+
+            //- Write
+            void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/kQRWallFunctions/kQRWallFunction/kQRWallFunctionFvPatchField.C

@@ -0,0 +1,155 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "kQRWallFunctionFvPatchField.H"
+#include "RASModel.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+#include "wallFvPatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+void kQRWallFunctionFvPatchField<Type>::checkType()
+{
+    if (!isA<wallFvPatch>(this->patch()))
+    {
+        FatalErrorIn("kQRWallFunctionFvPatchField::checkType()")
+            << "Invalid wall function specification" << nl
+            << "    Patch type for patch " << this->patch().name()
+            << " must be wall" << nl
+            << "    Current patch type is " << this->patch().type()
+            << nl << endl << abort(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+kQRWallFunctionFvPatchField<Type>::kQRWallFunctionFvPatchField
+(
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF
+)
+:
+    zeroGradientFvPatchField<Type>(p, iF)
+{
+    checkType();
+}
+
+
+template<class Type>
+kQRWallFunctionFvPatchField<Type>::kQRWallFunctionFvPatchField
+(
+    const kQRWallFunctionFvPatchField& ptf,
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    zeroGradientFvPatchField<Type>(ptf, p, iF, mapper)
+{
+    checkType();
+}
+
+
+template<class Type>
+kQRWallFunctionFvPatchField<Type>::kQRWallFunctionFvPatchField
+(
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    zeroGradientFvPatchField<Type>(p, iF, dict)
+{
+    checkType();
+}
+
+
+template<class Type>
+kQRWallFunctionFvPatchField<Type>::kQRWallFunctionFvPatchField
+(
+    const kQRWallFunctionFvPatchField& tkqrwfpf
+)
+:
+    zeroGradientFvPatchField<Type>(tkqrwfpf)
+{
+    checkType();
+}
+
+
+template<class Type>
+kQRWallFunctionFvPatchField<Type>::kQRWallFunctionFvPatchField
+(
+    const kQRWallFunctionFvPatchField& tkqrwfpf,
+    const DimensionedField<Type, volMesh>& iF
+)
+:
+    zeroGradientFvPatchField<Type>(tkqrwfpf, iF)
+{
+    checkType();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void kQRWallFunctionFvPatchField<Type>::evaluate
+(
+    const Pstream::commsTypes commsType
+)
+{
+    zeroGradientFvPatchField<Type>::evaluate(commsType);
+}
+
+
+template<class Type>
+void kQRWallFunctionFvPatchField<Type>::write(Ostream& os) const
+{
+    zeroGradientFvPatchField<Type>::write(os);
+    this->writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// ************************************************************************* //

src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/kQRWallFunctions/kQRWallFunction/kQRWallFunctionFvPatchField.H

@@ -0,0 +1,170 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::compressible::RASModels::kQRWallFunctionFvPatchField
+
+Description
+    Boundary condition for turbulence k, Q, and R when using wall functions.
+    Simply acts as a zero gradient condition.
+
+SourceFiles
+    kQRWallFunctionFvPatchField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef kQRWallFunctionFvPatchField_H
+#define kQRWallFunctionFvPatchField_H
+
+#include "zeroGradientFvPatchField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+                Class kQRWallFunctionFvPatchField Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class kQRWallFunctionFvPatchField
+:
+    public zeroGradientFvPatchField<Type>
+{
+
+    // Private member functions
+
+        //- Check the type of the patch
+        void checkType();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("kQRWallFunction");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        kQRWallFunctionFvPatchField
+        (
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        kQRWallFunctionFvPatchField
+        (
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        // kQRWallFunctionFvPatchField
+        //  onto a new patch
+        kQRWallFunctionFvPatchField
+        (
+            const kQRWallFunctionFvPatchField&,
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        kQRWallFunctionFvPatchField
+        (
+            const kQRWallFunctionFvPatchField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchField<Type> > clone() const
+        {
+            return tmp<fvPatchField<Type> >
+            (
+                new kQRWallFunctionFvPatchField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        kQRWallFunctionFvPatchField
+        (
+            const kQRWallFunctionFvPatchField&,
+            const DimensionedField<Type, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchField<Type> > clone
+        (
+            const DimensionedField<Type, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchField<Type> >
+            (
+                new kQRWallFunctionFvPatchField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Evaluation functions
+
+            //- Evaluate the patchField
+            virtual void evaluate
+            (
+                const Pstream::commsTypes commsType=Pstream::Pstream::blocking
+            );
+
+
+        // I-O
+
+            //- Write
+            void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "kQRWallFunctionFvPatchField.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/kQRWallFunctions/kQRWallFunction/kQRWallFunctionFvPatchFields.C

@@ -22,30 +22,30 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Global
-    wallOmega
-
-Description
-    Set wall dissipation in the omega matrix
-
 \*---------------------------------------------------------------------------*/
 
+#include "kQRWallFunctionFvPatchFields.H"
+#include "fvPatchFields.H"
+#include "addToRunTimeSelectionTable.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
 {
-    const fvPatchList& patches = mesh_.boundary();
 
-    forAll(patches, patchi)
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-    {
-        const fvPatch& p = patches[patchi];
 
-        if (isType<wallFvPatch>(p))
+makePatchFields(kQRWallFunction);
-        {
+
-            omegaEqn().setValues
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-            (
+
-                p.faceCells(),
+} // End namespace RASModels
-                omega_.boundaryField()[patchi].patchInternalField()
+} // End namespace compressible
-            );
+} // End namespace Foam
-        }
-    }
-}
 
 // ************************************************************************* //

src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/kQRWallFunctions/kQRWallFunction/kQRWallFunctionFvPatchFields.H

@@ -0,0 +1,56 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef kQRWallFunctionFvPatchFields_H
+#define kQRWallFunctionFvPatchFields_H
+
+#include "kQRWallFunctionFvPatchField.H"
+#include "fieldTypes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeFieldTypedefs(kQRWallFunction)
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutRoughWallFunction/mutRoughWallFunctionFvPatchScalarField.C

@@ -0,0 +1,248 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "mutRoughWallFunctionFvPatchScalarField.H"
+#include "RASModel.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+scalar mutRoughWallFunctionFvPatchScalarField::fnRough
+(
+    const scalar KsPlus,
+    const scalar Cs,
+    const scalar kappa
+) const
+{
+    // Set deltaB based on non-dimensional roughness height
+    scalar deltaB = 0.0;
+    if (KsPlus < 90.0)
+    {
+        deltaB =
+            1.0/kappa
+            *log((KsPlus - 2.25)/87.75 + Cs*KsPlus)
+            *sin(0.4258*(log(KsPlus) - 0.811));
+    }
+    else
+    {
+        deltaB = 1.0/kappa*log(1.0 + Cs*KsPlus);
+    }
+
+    return exp(min(deltaB*kappa, 50.0));
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+mutRoughWallFunctionFvPatchScalarField::
+mutRoughWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(p, iF),
+    Ks_(p.size(), 0.0),
+    Cs_(p.size(), 0.0)
+{}
+
+
+mutRoughWallFunctionFvPatchScalarField::
+mutRoughWallFunctionFvPatchScalarField
+(
+    const mutRoughWallFunctionFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
+    Ks_(ptf.Ks_, mapper),
+    Cs_(ptf.Cs_, mapper)
+{}
+
+
+mutRoughWallFunctionFvPatchScalarField::
+mutRoughWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchScalarField(p, iF, dict),
+    Ks_("Ks", dict, p.size()),
+    Cs_("Cs", dict, p.size())
+{}
+
+
+mutRoughWallFunctionFvPatchScalarField::
+mutRoughWallFunctionFvPatchScalarField
+(
+    const mutRoughWallFunctionFvPatchScalarField& nrwfpsf
+)
+:
+    fixedValueFvPatchScalarField(nrwfpsf),
+    Ks_(nrwfpsf.Ks_),
+    Cs_(nrwfpsf.Cs_)
+{}
+
+
+mutRoughWallFunctionFvPatchScalarField::
+mutRoughWallFunctionFvPatchScalarField
+(
+    const mutRoughWallFunctionFvPatchScalarField& nrwfpsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(nrwfpsf, iF),
+    Ks_(nrwfpsf.Ks_),
+    Cs_(nrwfpsf.Cs_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void mutRoughWallFunctionFvPatchScalarField::autoMap
+(
+    const fvPatchFieldMapper& m
+)
+{
+    fixedValueFvPatchScalarField::autoMap(m);
+    Ks_.autoMap(m);
+    Cs_.autoMap(m);
+}
+
+
+void mutRoughWallFunctionFvPatchScalarField::rmap
+(
+    const fvPatchScalarField& ptf,
+    const labelList& addr
+)
+{
+    fixedValueFvPatchScalarField::rmap(ptf, addr);
+
+    const mutRoughWallFunctionFvPatchScalarField& nrwfpsf =
+        refCast<const mutRoughWallFunctionFvPatchScalarField>(ptf);
+
+    Cs_.rmap(nrwfpsf.Cs_, addr);
+    Ks_.rmap(nrwfpsf.Ks_, addr);
+}
+
+
+void mutRoughWallFunctionFvPatchScalarField::updateCoeffs()
+{
+    const RASModel& ras = db().lookupObject<RASModel>("RASProperties");
+
+    const scalar Cmu = ras.Cmu().value();
+    const scalar Cmu25 = pow(Cmu, 0.25);
+    const scalar kappa = ras.kappa().value();
+    const scalar E = ras.E().value();
+    scalar yPlusLam = ras.yPlusLam();
+
+    const scalarField& y = ras.y()[patch().index()];
+
+    const scalarField& rhow =
+        patch().lookupPatchField<volScalarField, scalar>("rho");
+
+    const scalarField& k = db().lookupObject<volScalarField>("k");
+
+    const scalarField& muw =
+        patch().lookupPatchField<volScalarField, scalar>("mu");
+
+    scalarField& mutw = *this;
+
+    forAll(mutw, faceI)
+    {
+        label faceCellI = patch().faceCells()[faceI];
+
+        scalar uStar = Cmu25*sqrt(k[faceCellI]);
+
+        scalar yPlus = uStar*y[faceI]/(muw[faceI]/rhow[faceI]);
+
+        scalar KsPlus = uStar*Ks_[faceI]/(muw[faceI]/rhow[faceI]);
+
+        scalar Edash = E;
+        scalar yPlusLamNew = yPlusLam;
+        if (KsPlus > 2.25)
+        {
+            Edash = E/fnRough(KsPlus, Cs_[faceI], kappa);
+            yPlusLam = ras.yPlusLam(kappa, Edash);
+        }
+
+        if (debug)
+        {
+            Info<< "yPlus = " << yPlus
+                << ", KsPlus = " << KsPlus
+                << ", Edash = " << Edash
+                << ", yPlusLam = " << yPlusLam
+                << endl;
+        }
+
+        if (yPlus > yPlusLamNew)
+        {
+            mutw[faceI] = muw[faceI]*(yPlus*kappa/log(Edash*yPlus) - 1);
+        }
+        else
+        {
+            mutw[faceI] = 0.0;
+        }
+    }
+}
+
+
+void mutRoughWallFunctionFvPatchScalarField::write(Ostream& os) const
+{
+    fvPatchField<scalar>::write(os);
+    Cs_.writeEntry("Cs", os);
+    Ks_.writeEntry("Ks", os);
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField(fvPatchScalarField, mutRoughWallFunctionFvPatchScalarField);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// ************************************************************************* //

src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutRoughWallFunction/mutRoughWallFunctionFvPatchScalarField.H

@@ -0,0 +1,194 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::compressible::RASModels::nutRoughWallFunctionFvPatchScalarField
+
+Description
+    Boundary condition for turbulent (kinematic) viscosity when using wall
+    functions for rough walls.
+
+    Manipulates the E parameter to account for roughness effects, based on
+    KsPlus.
+
+    - roughness height = sand-grain roughness (0 for smooth walls)
+    - roughness constant = 0.5-1.0 (0.5 default)
+
+SourceFiles
+    mutRoughWallFunctionFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef mutRoughWallFunctionFvPatchScalarField_H
+#define mutRoughWallFunctionFvPatchScalarField_H
+
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+          Class mutRoughWallFunctionFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class mutRoughWallFunctionFvPatchScalarField
+:
+    public fixedValueFvPatchScalarField
+{
+    // Private data
+
+        //- Roughness height
+        scalarField Ks_;
+
+        //- Roughness constant
+        scalarField Cs_;
+
+
+    // Private member functions
+
+        //- Compute the roughness function
+        scalar fnRough
+        (
+            const scalar KsPlus,
+            const scalar Cs,
+            const scalar kappa
+        ) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("mutRoughWallFunction");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        mutRoughWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        mutRoughWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  mutRoughWallFunctionFvPatchScalarField
+        //  onto a new patch
+        mutRoughWallFunctionFvPatchScalarField
+        (
+            const mutRoughWallFunctionFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        mutRoughWallFunctionFvPatchScalarField
+        (
+            const mutRoughWallFunctionFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new mutRoughWallFunctionFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        mutRoughWallFunctionFvPatchScalarField
+        (
+            const mutRoughWallFunctionFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new mutRoughWallFunctionFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Mapping functions
+
+            //- Map (and resize as needed) from self given a mapping object
+            virtual void autoMap
+            (
+                const fvPatchFieldMapper&
+            );
+
+            //- Reverse map the given fvPatchField onto this fvPatchField
+            virtual void rmap
+            (
+                const fvPatchScalarField&,
+                const labelList&
+            );
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+
+        // I-O
+
+            //- Write
+            void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardRoughWallFunction/mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C

@@ -24,7 +24,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "mutStandardRoughWallFunctionFvPatchScalarField.H"
+#include "mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.H"
 #include "RASModel.H"
 #include "fvPatchFieldMapper.H"
 #include "volFields.H"
@@ -41,8 +41,8 @@ namespace RASModels
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-mutStandardRoughWallFunctionFvPatchScalarField::
+mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::
-mutStandardRoughWallFunctionFvPatchScalarField
+mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
 (
     const fvPatch& p,
     const DimensionedField<scalar, volMesh>& iF
@@ -55,10 +55,10 @@ mutStandardRoughWallFunctionFvPatchScalarField
 {}
 
 
-mutStandardRoughWallFunctionFvPatchScalarField::
+mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::
-mutStandardRoughWallFunctionFvPatchScalarField
+mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
 (
-    const mutStandardRoughWallFunctionFvPatchScalarField& ptf,
+    const mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField& ptf,
     const fvPatch& p,
     const DimensionedField<scalar, volMesh>& iF,
     const fvPatchFieldMapper& mapper
@@ -71,8 +71,8 @@ mutStandardRoughWallFunctionFvPatchScalarField
 {}
 
 
-mutStandardRoughWallFunctionFvPatchScalarField::
+mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::
-mutStandardRoughWallFunctionFvPatchScalarField
+mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
 (
     const fvPatch& p,
     const DimensionedField<scalar, volMesh>& iF,
@@ -86,10 +86,10 @@ mutStandardRoughWallFunctionFvPatchScalarField
 {}
 
 
-mutStandardRoughWallFunctionFvPatchScalarField::
+mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::
-mutStandardRoughWallFunctionFvPatchScalarField
+mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
 (
-    const mutStandardRoughWallFunctionFvPatchScalarField& tppsf
+    const mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField& tppsf
 )
 :
     fixedValueFvPatchScalarField(tppsf),
@@ -99,10 +99,10 @@ mutStandardRoughWallFunctionFvPatchScalarField
 {}
 
 
-mutStandardRoughWallFunctionFvPatchScalarField::
+mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::
-mutStandardRoughWallFunctionFvPatchScalarField
+mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
 (
-    const mutStandardRoughWallFunctionFvPatchScalarField& tppsf,
+    const mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField& tppsf,
     const DimensionedField<scalar, volMesh>& iF
 )
 :
@@ -115,7 +115,7 @@ mutStandardRoughWallFunctionFvPatchScalarField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void mutStandardRoughWallFunctionFvPatchScalarField::evaluate
+void mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::evaluate
 (
     const Pstream::commsTypes
 )
@@ -282,7 +282,10 @@ void mutStandardRoughWallFunctionFvPatchScalarField::evaluate
 }
 
 
-void mutStandardRoughWallFunctionFvPatchScalarField::write(Ostream& os) const
+void mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::write
+(
+    Ostream& os
+) const
 {
     fixedValueFvPatchScalarField::write(os);
     os.writeKeyword("roughnessHeight")
@@ -299,7 +302,7 @@ void mutStandardRoughWallFunctionFvPatchScalarField::write(Ostream& os) const
 makePatchTypeField
 (
     fvPatchScalarField,
-    mutStandardRoughWallFunctionFvPatchScalarField
+    mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
 );
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardRoughWallFunction/mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.H

@@ -23,18 +23,19 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::compressible::RASModels::mutStandardRoughWallFunctionFvPatchScalarField
+    Foam::compressible::RASModels::
+        mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
 
 Description
     Wall function boundary condition for rough walls
 
 SourceFiles
-    mutStandardWallFunctionFvPatchScalarField.C
+    mutSpalartAllamarasStandardWallFunctionFvPatchScalarField.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef mutStandardRoughWallFunctionFvPatchScalarField_H
+#ifndef mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField_H
-#define mutStandardRoughWallFunctionFvPatchScalarField_H
+#define mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField_H
 
 #include "fixedValueFvPatchFields.H"
 
@@ -48,10 +49,10 @@ namespace RASModels
 {
 
 /*---------------------------------------------------------------------------*\
-                    Class mutWallFunctionFvPatch Declaration
+Class mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField Declaration
 \*---------------------------------------------------------------------------*/
 
-class mutStandardRoughWallFunctionFvPatchScalarField
+class mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
 :
     public fixedValueFvPatchScalarField
 {
@@ -65,41 +66,41 @@ class mutStandardRoughWallFunctionFvPatchScalarField
 public:
 
     //- Runtime type information
-    TypeName("mutStandardRoughWallFunction");
+    TypeName("mutSpalartAllmarasStandardRoughWallFunction");
 
 
     // Constructors
 
         //- Construct from patch and internal field
-        mutStandardRoughWallFunctionFvPatchScalarField
+        mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
         (
             const fvPatch&,
             const DimensionedField<scalar, volMesh>&
         );
 
         //- Construct from patch, internal field and dictionary
-        mutStandardRoughWallFunctionFvPatchScalarField
+        mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
         (
             const fvPatch&,
             const DimensionedField<scalar, volMesh>&,
             const dictionary&
         );
 
-        //- Construct by mapping given 
+        //- Construct by mapping given
-        //  mutStandardRoughWallFunctionFvPatchScalarField
+        //  mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
         //  onto a new patch
-        mutStandardRoughWallFunctionFvPatchScalarField
+        mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
         (
-            const mutStandardRoughWallFunctionFvPatchScalarField&,
+            const mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField&,
             const fvPatch&,
             const DimensionedField<scalar, volMesh>&,
             const fvPatchFieldMapper&
         );
 
         //- Construct as copy
-        mutStandardRoughWallFunctionFvPatchScalarField
+        mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
         (
-            const mutStandardRoughWallFunctionFvPatchScalarField&
+            const mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField&
         );
 
         //- Construct and return a clone
@@ -107,14 +108,17 @@ public:
         {
             return tmp<fvPatchScalarField>
             (
-                new mutStandardRoughWallFunctionFvPatchScalarField(*this)
+                new mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
+                (
+                    *this
+                )
             );
         }
 
         //- Construct as copy setting internal field reference
-        mutStandardRoughWallFunctionFvPatchScalarField
+        mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
         (
-            const mutStandardRoughWallFunctionFvPatchScalarField&,
+            const mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField&,
             const DimensionedField<scalar, volMesh>&
         );
 
@@ -126,7 +130,11 @@ public:
         {
             return tmp<fvPatchScalarField>
             (
-                new mutStandardRoughWallFunctionFvPatchScalarField(*this, iF)
+                new mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
+                (
+                   *this,
+                    iF
+                )
             );
         }
 

src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardWallFunction/mutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C

@@ -0,0 +1,175 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "mutSpalartAllmarasStandardWallFunctionFvPatchScalarField.H"
+#include "RASModel.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::
+mutSpalartAllmarasStandardWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(p, iF)
+{}
+
+
+mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::
+mutSpalartAllmarasStandardWallFunctionFvPatchScalarField
+(
+    const mutSpalartAllmarasStandardWallFunctionFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchScalarField(ptf, p, iF, mapper)
+{}
+
+
+mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::
+mutSpalartAllmarasStandardWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchScalarField(p, iF, dict)
+{}
+
+
+mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::
+mutSpalartAllmarasStandardWallFunctionFvPatchScalarField
+(
+    const mutSpalartAllmarasStandardWallFunctionFvPatchScalarField& tppsf
+)
+:
+    fixedValueFvPatchScalarField(tppsf)
+{}
+
+
+mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::
+mutSpalartAllmarasStandardWallFunctionFvPatchScalarField
+(
+    const mutSpalartAllmarasStandardWallFunctionFvPatchScalarField& tppsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(tppsf, iF)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::evaluate
+(
+    const Pstream::commsTypes
+)
+{
+    const RASModel& ras = db().lookupObject<RASModel>("RASProperties");
+
+    scalar kappa = ras.kappa().value();
+    scalar E = ras.E().value();
+    scalar yPlusLam = ras.yPlusLam();
+
+    const scalarField& ry = patch().deltaCoeffs();
+
+    const fvPatchVectorField& U =
+        patch().lookupPatchField<volVectorField, vector>("U");
+
+    scalarField magUp = mag(U.patchInternalField() - U);
+
+    const scalarField& rhow =
+        patch().lookupPatchField<volScalarField, scalar>("rho");
+
+    const scalarField& muw =
+        patch().lookupPatchField<volScalarField, scalar>("mu");
+    scalarField& mutw = *this;
+
+    scalarField magFaceGradU = mag(U.snGrad());
+
+    forAll(mutw, faceI)
+    {
+        scalar magUpara = magUp[faceI];
+
+        scalar kappaRe = kappa*magUpara/((muw[faceI]/rhow[faceI])*ry[faceI]);
+
+        scalar yPlus = yPlusLam;
+        scalar ryPlusLam = 1.0/yPlus;
+
+        int iter = 0;
+        scalar yPlusLast = 0.0;
+
+        do
+        {
+            yPlusLast = yPlus;
+            yPlus = (kappaRe + yPlus)/(1.0 + log(E*yPlus));
+
+        } while(mag(ryPlusLam*(yPlus - yPlusLast)) > 0.01 && ++iter < 10 );
+
+        if (yPlus > yPlusLam)
+        {
+            mutw[faceI] = muw[faceI]*(yPlus*kappa/log(E*yPlus) - 1);
+        }
+        else
+        {
+            mutw[faceI] = 0.0;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField
+(
+    fvPatchScalarField,
+    mutSpalartAllmarasStandardWallFunctionFvPatchScalarField
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// ************************************************************************* //

src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardWallFunction/mutSpalartAllmarasStandardWallFunctionFvPatchScalarField.H

@@ -0,0 +1,158 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::compressible::RASModels::
+    mutSpalartAllmarasStandardWallFunctionFvPatchScalarField
+
+Description
+    Wall function boundary condition for walls
+
+SourceFiles
+    mutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef mutSpalartAllmarasStandardWallFunctionFvPatchScalarField_H
+#define mutSpalartAllmarasStandardWallFunctionFvPatchScalarField_H
+
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+  Class mutSpalartAllmarasStandardWallFunctionFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class mutSpalartAllmarasStandardWallFunctionFvPatchScalarField
+:
+    public fixedValueFvPatchScalarField
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("mutSpalartAllmarasStandardWallFunction");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        mutSpalartAllmarasStandardWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        mutSpalartAllmarasStandardWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  mutSpalartAllmarasStandardWallFunctionFvPatchScalarField
+        //  onto a new patch
+        mutSpalartAllmarasStandardWallFunctionFvPatchScalarField
+        (
+            const mutSpalartAllmarasStandardWallFunctionFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        mutSpalartAllmarasStandardWallFunctionFvPatchScalarField
+        (
+            const mutSpalartAllmarasStandardWallFunctionFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new mutSpalartAllmarasStandardWallFunctionFvPatchScalarField
+                (
+                    *this
+                )
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        mutSpalartAllmarasStandardWallFunctionFvPatchScalarField
+        (
+            const mutSpalartAllmarasStandardWallFunctionFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new mutSpalartAllmarasStandardWallFunctionFvPatchScalarField
+                (
+                    *this,
+                    iF
+                )
+            );
+        }
+
+
+    // Member functions
+
+        // Evaluation functions
+
+            //- Evaluate the patchField
+            virtual void evaluate
+            (
+                const Pstream::commsTypes commsType=Pstream::blocking
+            );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasWallFunction/mutSpalartAllmarasWallFunctionFvPatchScalarField.C

@@ -0,0 +1,187 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "mutSpalartAllmarasWallFunctionFvPatchScalarField.H"
+#include "RASModel.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+mutSpalartAllmarasWallFunctionFvPatchScalarField::
+mutSpalartAllmarasWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(p, iF)
+{}
+
+
+mutSpalartAllmarasWallFunctionFvPatchScalarField::
+mutSpalartAllmarasWallFunctionFvPatchScalarField
+(
+    const mutSpalartAllmarasWallFunctionFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchScalarField(ptf, p, iF, mapper)
+{}
+
+
+mutSpalartAllmarasWallFunctionFvPatchScalarField::
+mutSpalartAllmarasWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchScalarField(p, iF, dict)
+{}
+
+
+mutSpalartAllmarasWallFunctionFvPatchScalarField::
+mutSpalartAllmarasWallFunctionFvPatchScalarField
+(
+    const mutSpalartAllmarasWallFunctionFvPatchScalarField& tppsf
+)
+:
+    fixedValueFvPatchScalarField(tppsf)
+{}
+
+
+mutSpalartAllmarasWallFunctionFvPatchScalarField::
+mutSpalartAllmarasWallFunctionFvPatchScalarField
+(
+    const mutSpalartAllmarasWallFunctionFvPatchScalarField& tppsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(tppsf, iF)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void mutSpalartAllmarasWallFunctionFvPatchScalarField::evaluate
+(
+    const Pstream::commsTypes
+)
+{
+    const RASModel& ras = db().lookupObject<RASModel>("RASProperties");
+
+    scalar kappa = ras.kappa().value();
+    scalar E = ras.E().value();
+
+    const scalarField& ry = patch().deltaCoeffs();
+
+    const fvPatchVectorField& U =
+        patch().lookupPatchField<volVectorField, vector>("U");
+
+    scalarField magUp = mag(U.patchInternalField() - U);
+
+    const scalarField& rhow =
+        patch().lookupPatchField<volScalarField, scalar>("rho");
+
+    const scalarField& muw =
+        patch().lookupPatchField<volScalarField, scalar>("mu");
+
+    scalarField& mutw = *this;
+
+    scalarField magFaceGradU = mag(U.snGrad());
+
+    forAll(mutw, faceI)
+    {
+        scalar magUpara = magUp[faceI];
+
+        scalar utau =
+            sqrt((mutw[faceI] + muw[faceI])*magFaceGradU[faceI]/rhow[faceI]);
+
+        if (utau > VSMALL)
+        {
+            int iter = 0;
+            scalar err = GREAT;
+
+            do
+            {
+                scalar kUu = min(kappa*magUpara/utau, 50);
+                scalar fkUu = exp(kUu) - 1 - kUu*(1 + 0.5*kUu);
+
+                scalar f =
+                    - utau/(ry[faceI]*(muw[faceI]/rhow[faceI]))
+                    + magUpara/utau
+                    + 1/E*(fkUu - 1.0/6.0*kUu*sqr(kUu));
+
+                scalar df =
+                    1.0/(ry[faceI]*(muw[faceI]/rhow[faceI]))
+                  + magUpara/sqr(utau)
+                  + 1/E*kUu*fkUu/utau;
+
+                scalar utauNew = utau + f/df;
+                err = mag((utau - utauNew)/utau);
+                utau = utauNew;
+
+            } while (utau > VSMALL && err > 0.01 && ++iter < 10);
+
+            mutw[faceI] = max
+            (
+                rhow[faceI]*sqr(max(utau, 0))/magFaceGradU[faceI]- muw[faceI],
+                0.0
+            );
+        }
+        else
+        {
+            mutw[faceI] = 0;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField(fvPatchScalarField, mutSpalartAllmarasWallFunctionFvPatchScalarField);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// ************************************************************************* //

src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasWallFunction/mutSpalartAllmarasWallFunctionFvPatchScalarField.H

@@ -0,0 +1,150 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::compressible::RASModels::
+    mutSpalartAllmarasWallFunctionFvPatchScalarField
+
+Description
+    Wall function boundary condition for walls
+
+SourceFiles
+    mutSpalartAllmarasWallFunctionFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef mutSpalartAllmarasWallFunctionFvPatchScalarField_H
+#define mutSpalartAllmarasWallFunctionFvPatchScalarField_H
+
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+      Class mutSpalartAllmarasWallFunctionFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class mutSpalartAllmarasWallFunctionFvPatchScalarField
+:
+    public fixedValueFvPatchScalarField
+{
+public:
+
+    //- Runtime type information
+    TypeName("mutSpalartAllmarasWallFunction");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        mutSpalartAllmarasWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        mutSpalartAllmarasWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  mutSpalartAllmarasWallFunctionFvPatchScalarField
+        //  onto a new patch
+        mutSpalartAllmarasWallFunctionFvPatchScalarField
+        (
+            const mutSpalartAllmarasWallFunctionFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        mutSpalartAllmarasWallFunctionFvPatchScalarField
+        (
+            const mutSpalartAllmarasWallFunctionFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new mutSpalartAllmarasWallFunctionFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        mutSpalartAllmarasWallFunctionFvPatchScalarField
+        (
+            const mutSpalartAllmarasWallFunctionFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new mutSpalartAllmarasWallFunctionFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Evaluation functions
+
+            //- Evaluate the patchField
+            virtual void evaluate
+            (
+                const Pstream::commsTypes commsType=Pstream::blocking
+            );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutWallFunction/mutWallFunctionFvPatchScalarField.C

@@ -0,0 +1,162 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "mutWallFunctionFvPatchScalarField.H"
+#include "RASModel.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+mutWallFunctionFvPatchScalarField::
+mutWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(p, iF)
+{}
+
+
+mutWallFunctionFvPatchScalarField::
+mutWallFunctionFvPatchScalarField
+(
+    const mutWallFunctionFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchScalarField(ptf, p, iF, mapper)
+{}
+
+
+mutWallFunctionFvPatchScalarField::
+mutWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchScalarField(p, iF, dict)
+{}
+
+
+mutWallFunctionFvPatchScalarField::
+mutWallFunctionFvPatchScalarField
+(
+    const mutWallFunctionFvPatchScalarField& tppsf
+)
+:
+    fixedValueFvPatchScalarField(tppsf)
+{}
+
+
+mutWallFunctionFvPatchScalarField::
+mutWallFunctionFvPatchScalarField
+(
+    const mutWallFunctionFvPatchScalarField& tppsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(tppsf, iF)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void mutWallFunctionFvPatchScalarField::updateCoeffs()
+{
+    const RASModel& ras = db().lookupObject<RASModel>("RASProperties");
+
+    const scalar Cmu = ras.Cmu().value();
+    const scalar Cmu25 = pow(Cmu, 0.25);
+    const scalar kappa = ras.kappa().value();
+    const scalar E = ras.E().value();
+    const scalar yPlusLam = ras.yPlusLam();
+
+    const scalarField& y = ras.y()[patch().index()];
+
+    const scalarField& rhow =
+        patch().lookupPatchField<volScalarField, scalar>("rho");
+
+    const volScalarField& k = db().lookupObject<volScalarField>("k");
+
+    const scalarField& muw =
+        patch().lookupPatchField<volScalarField, scalar>("mu");
+
+    scalarField& mutw = *this;
+
+    forAll(mutw, faceI)
+    {
+        label faceCellI = patch().faceCells()[faceI];
+
+        scalar yPlus =
+            Cmu25*y[faceI]*sqrt(k[faceCellI])
+           /(muw[faceI]/rhow[faceI]);
+
+        if (yPlus > yPlusLam)
+        {
+            mutw[faceI] = muw[faceI]*(yPlus*kappa/log(E*yPlus) - 1);
+        }
+        else
+        {
+            mutw[faceI] = 0.0;
+        }
+    }
+}
+
+
+void mutWallFunctionFvPatchScalarField::write(Ostream& os) const
+{
+    fvPatchField<scalar>::write(os);
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField(fvPatchScalarField, mutWallFunctionFvPatchScalarField);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// ************************************************************************* //

src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutWallFunction/mutWallFunctionFvPatchScalarField.H

@@ -0,0 +1,155 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::compressible::RASModels::mutWallFunctionFvPatchScalarField
+
+Description
+    Boundary condition for turbulent (kinematic) viscosity when using wall
+    functions
+    - replicates OpenFOAM v1.5 (and earlier) behaviour
+
+SourceFiles
+    mutWallFunctionFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef mutWallFunctionFvPatchScalarField_H
+#define mutWallFunctionFvPatchScalarField_H
+
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+             Class mutWallFunctionFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class mutWallFunctionFvPatchScalarField
+:
+    public fixedValueFvPatchScalarField
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("mutWallFunction");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        mutWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        mutWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  mutWallFunctionFvPatchScalarField
+        //  onto a new patch
+        mutWallFunctionFvPatchScalarField
+        (
+            const mutWallFunctionFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        mutWallFunctionFvPatchScalarField
+        (
+            const mutWallFunctionFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new mutWallFunctionFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        mutWallFunctionFvPatchScalarField
+        (
+            const mutWallFunctionFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new mutWallFunctionFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+
+        // I-O
+
+            //- Write
+            void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C

@@ -0,0 +1,209 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "omegaWallFunctionFvPatchScalarField.H"
+#include "RASModel.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+#include "wallFvPatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void omegaWallFunctionFvPatchScalarField::checkType()
+{
+    if (!isA<wallFvPatch>(patch()))
+    {
+        FatalErrorIn("omegaWallFunctionFvPatchScalarField::checkType()")
+            << "Invalid wall function specification" << nl
+            << "    Patch type for patch " << patch().name()
+            << " must be wall" << nl
+            << "    Current patch type is " << patch().type() << nl << endl
+            << abort(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedInternalValueFvPatchField<scalar>(p, iF)
+{
+    checkType();
+}
+
+
+omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
+(
+    const omegaWallFunctionFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedInternalValueFvPatchField<scalar>(ptf, p, iF, mapper)
+{
+    checkType();
+}
+
+
+omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedInternalValueFvPatchField<scalar>(p, iF, dict)
+{
+    checkType();
+}
+
+
+omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
+(
+    const omegaWallFunctionFvPatchScalarField& ewfpsf
+)
+:
+    fixedInternalValueFvPatchField<scalar>(ewfpsf)
+{
+    checkType();
+}
+
+
+omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
+(
+    const omegaWallFunctionFvPatchScalarField& ewfpsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedInternalValueFvPatchField<scalar>(ewfpsf, iF)
+{
+    checkType();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void omegaWallFunctionFvPatchScalarField::updateCoeffs()
+{
+    const RASModel& ras = db().lookupObject<RASModel>("RASProperties");
+
+    const scalar Cmu = ras.Cmu().value();
+    const scalar Cmu25 = pow(Cmu, 0.25);
+
+    const scalar kappa = ras.kappa().value();
+    const scalar yPlusLam = ras.yPlusLam();
+
+    const scalarField& y = ras.y()[patch().index()];
+
+    volScalarField& G = const_cast<volScalarField&>
+        (db().lookupObject<volScalarField>("G"));
+
+    volScalarField& omega = const_cast<volScalarField&>
+        (db().lookupObject<volScalarField>("omega"));
+
+    const scalarField& k = db().lookupObject<volScalarField>("k");
+
+    const scalarField& rhow =
+        patch().lookupPatchField<volScalarField, scalar>("rho");
+
+    const scalarField& muw =
+        patch().lookupPatchField<volScalarField, scalar>("mu");
+
+    const scalarField& mutw =
+        patch().lookupPatchField<volScalarField, scalar>("mut");
+
+    const fvPatchVectorField& Uw =
+        patch().lookupPatchField<volVectorField, vector>("U");
+
+    const scalarField magGradUw = mag(Uw.snGrad());
+
+    // Set epsilon and G
+    forAll(mutw, faceI)
+    {
+        label faceCellI = patch().faceCells()[faceI];
+
+        scalar yPlus =
+            Cmu25*y[faceI]*sqrt(k[faceCellI])
+           /(muw[faceI]/rhow[faceI]);
+
+        omega[faceCellI] = sqrt(k[faceCellI])/(Cmu25*kappa*y[faceI]);
+
+        if (yPlus > yPlusLam)
+        {
+            G[faceCellI] =
+                (mutw[faceI] + muw[faceI])
+               *magGradUw[faceI]
+               *Cmu25*sqrt(k[faceCellI])
+               /(kappa*y[faceI]);
+        }
+        else
+        {
+            G[faceCellI] = 0.0;
+        }
+    }
+
+    // TODO: perform averaging for cells sharing more than one boundary face
+}
+
+
+void omegaWallFunctionFvPatchScalarField::write(Ostream& os) const
+{
+    fixedInternalValueFvPatchField<scalar>::write(os);
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField
+(
+    fvPatchScalarField,
+    omegaWallFunctionFvPatchScalarField
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// ************************************************************************* //

src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.H

@@ -0,0 +1,159 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::compressible::RASModels::omegaWallFunctionFvPatchScalarField
+
+Description
+    Replaces functionality in wallFunctionsI.H
+
+SourceFiles
+    omegaWallFunctionFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef omegaWallFunctionFvPatchScalarField_H
+#define omegaWallFunctionFvPatchScalarField_H
+
+#include "fixedInternalValueFvPatchField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+           Class omegaWallFunctionFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class omegaWallFunctionFvPatchScalarField
+:
+    public fixedInternalValueFvPatchField<scalar>
+{
+
+    // Private member functions
+
+        //- Check the type of the patch
+        void checkType();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("omegaWallFunction");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        omegaWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        omegaWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        // omegaWallFunctionFvPatchScalarField
+        //  onto a new patch
+        omegaWallFunctionFvPatchScalarField
+        (
+            const omegaWallFunctionFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        omegaWallFunctionFvPatchScalarField
+        (
+            const omegaWallFunctionFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new omegaWallFunctionFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        omegaWallFunctionFvPatchScalarField
+        (
+            const omegaWallFunctionFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new omegaWallFunctionFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+
+        // I-O
+
+            //- Write
+            void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/turbulenceModels/RAS/compressible/kEpsilon/kEpsilon.C

@@ -28,6 +28,8 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "wallFvPatch.H"
 
+#include "backwardsCompatibilityWallFunctions.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -125,12 +127,11 @@ kEpsilon::kEpsilon
             "k",
             runTime_.timeName(),
             mesh_,
-            IOobject::MUST_READ,
+            IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh_
+        autoCreateK("k", mesh_)
     ),
-
     epsilon_
     (
         IOobject
@@ -138,12 +139,11 @@ kEpsilon::kEpsilon
             "epsilon",
             runTime_.timeName(),
             mesh_,
-            IOobject::MUST_READ,
+            IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh_
+        autoCreateEpsilon("epsilon", mesh_)
     ),
-
     mut_
     (
         IOobject
@@ -152,12 +152,13 @@ kEpsilon::kEpsilon
             runTime_.timeName(),
             mesh_,
             IOobject::NO_READ,
-            IOobject::NO_WRITE
+            IOobject::AUTO_WRITE
         ),
-        Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_)
+        autoCreateMut("mut", mesh_)
     )
 {
-#   include "wallViscosityI.H"
+    mut_ = Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    mut_.correctBoundaryConditions();
 
     printCoeffs();
 }
@@ -243,7 +244,7 @@ void kEpsilon::correct()
     {
         // Re-calculate viscosity
         mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
-#       include "wallViscosityI.H"
+        mut_.correctBoundaryConditions();
         return;
     }
 
@@ -257,10 +258,11 @@ void kEpsilon::correct()
     }
 
     tmp<volTensorField> tgradU = fvc::grad(U_);
-    volScalarField G = mut_*(tgradU() && dev(twoSymm(tgradU())));
+    volScalarField G("G", mut_*(tgradU() && dev(twoSymm(tgradU()))));
     tgradU.clear();
 
-#   include "wallFunctionsI.H"
+    // Update espsilon and G at the wall
+    epsilon_.boundaryField().updateCoeffs();
 
     // Dissipation equation
     tmp<fvScalarMatrix> epsEqn
@@ -274,10 +276,10 @@ void kEpsilon::correct()
       - fvm::Sp(C2_*rho_*epsilon_/k_, epsilon_)
     );
 
-#   include "wallDissipationI.H"
-
     epsEqn().relax();
 
+    epsEqn().boundaryManipulate(epsilon_.boundaryField());
+
     solve(epsEqn);
     bound(epsilon_, epsilon0_);
 
@@ -302,9 +304,7 @@ void kEpsilon::correct()
 
     // Re-calculate viscosity
     mut_ = rho_*Cmu_*sqr(k_)/epsilon_;
-
+    mut_.correctBoundaryConditions();
-#   include "wallViscosityI.H"
-
 }
 
 

src/turbulenceModels/RAS/compressible/kOmegaSST/checkkOmegaPatchFieldTypes.H

@@ -1,19 +0,0 @@
-if (!isType<zeroGradientFvPatchScalarField>(k_.boundaryField()[patchi]))
-{
-    FatalErrorIn("wall-function evaluation")
-        << k_.boundaryField()[patchi].type()
-        << " is the wrong k patchField type for wall-functions on patch "
-        << curPatch.name() << nl
-        << "    should be zeroGradient"
-        << exit(FatalError);
-}
-
-if (!isType<zeroGradientFvPatchScalarField>(omega_.boundaryField()[patchi]))
-{
-    FatalErrorIn("wall-function evaluation")
-        << omega_.boundaryField()[patchi].type()
-        << " is the wrong omega patchField type for wall-functions on patch "
-        << curPatch.name() << nl
-        << "    should be zeroGradient"
-        << exit(FatalError);
-}

src/turbulenceModels/RAS/compressible/kOmegaSST/kOmegaSST.C

@@ -28,6 +28,8 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "wallFvPatch.H"
 
+#include "backwardsCompatibilityWallFunctions.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -228,12 +230,11 @@ kOmegaSST::kOmegaSST
             "k",
             runTime_.timeName(),
             mesh_,
-            IOobject::MUST_READ,
+            IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh_
+        autoCreateK("k", mesh_)
     ),
-
     omega_
     (
         IOobject
@@ -241,12 +242,11 @@ kOmegaSST::kOmegaSST
             "omega",
             runTime_.timeName(),
             mesh_,
-            IOobject::MUST_READ,
+            IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh_
+        autoCreateOmega("omega", mesh_)
     ),
-
     mut_
     (
         IOobject
@@ -255,12 +255,13 @@ kOmegaSST::kOmegaSST
             runTime_.timeName(),
             mesh_,
             IOobject::NO_READ,
-            IOobject::NO_WRITE
+            IOobject::AUTO_WRITE
         ),
-        a1_*rho_*k_/max(a1_*omega_, F2()*sqrt(magSqr(symm(fvc::grad(U_)))))
+        autoCreateMut("mut", mesh_)
     )
 {
-#   include "kOmegaWallViscosityI.H"
+    mut_ = a1_*rho_*k_/max(a1_*omega_, F2()*sqrt(magSqr(symm(fvc::grad(U_)))));
+    mut_.correctBoundaryConditions();
 
     printCoeffs();
 }
@@ -353,7 +354,8 @@ void kOmegaSST::correct()
         mut_ =
             a1_*rho_*k_
            /max(a1_*omega_, F2()*sqrt(magSqr(symm(fvc::grad(U_)))));
-#       include "kOmegaWallViscosityI.H"
+        mut_.correctBoundaryConditions();
+
         return;
     }
 
@@ -374,10 +376,11 @@ void kOmegaSST::correct()
     tmp<volTensorField> tgradU = fvc::grad(U_);
     volScalarField S2 = magSqr(symm(tgradU()));
     volScalarField GbyMu = (tgradU() && dev(twoSymm(tgradU())));
-    volScalarField G = mut_*GbyMu;
+    volScalarField G("G", mut_*GbyMu);
     tgradU.clear();
 
-#   include "kOmegaWallFunctionsI.H"
+    // Update omega and G at the wall
+    omega_.boundaryField().updateCoeffs();
 
     volScalarField CDkOmega =
         (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_;
@@ -404,7 +407,7 @@ void kOmegaSST::correct()
 
     omegaEqn().relax();
 
-#   include "wallOmegaI.H"
+    omegaEqn().boundaryManipulate(omega_.boundaryField());
 
     solve(omegaEqn);
     bound(omega_, omega0_);
@@ -428,9 +431,7 @@ void kOmegaSST::correct()
 
     // Re-calculate viscosity
     mut_ = a1_*rho_*k_/max(a1_*omega_, F2()*sqrt(S2));
-
+    mut_.correctBoundaryConditions();
-#   include "kOmegaWallViscosityI.H"
-
 }
 
 

src/turbulenceModels/RAS/compressible/kOmegaSST/kOmegaWallFunctionsI.H

@@ -1,127 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Global
-    kOmegaWallFunctions
-
-Description
-    Calculate wall generation and frequency omega from wall-functions.
-
-\*---------------------------------------------------------------------------*/
-
-{
-    labelList cellBoundaryFaceCount(omega_.size(), 0);
-
-    scalar Cmu25 = pow(Cmu_.value(), 0.25);
-
-    const fvPatchList& patches = mesh_.boundary();
-
-    //- Initialise the near-wall omega and G fields to zero
-    forAll(patches, patchi)
-    {
-        const fvPatch& curPatch = patches[patchi];
-
-        if (isType<wallFvPatch>(curPatch))
-        {
-            forAll(curPatch, facei)
-            {
-                label faceCelli = curPatch.faceCells()[facei];
-
-                omega_[faceCelli] = 0.0;
-                G[faceCelli] = 0.0;
-            }
-        }
-    }
-
-    //- Accumulate the wall face contributions to omega and G
-    //  Increment cellBoundaryFaceCount for each face for averaging
-    forAll(patches, patchi)
-    {
-        const fvPatch& curPatch = patches[patchi];
-
-        if (isType<wallFvPatch>(curPatch))
-        {
-#           include "checkkOmegaPatchFieldTypes.H"
-
-            const scalarField& rhow = rho_.boundaryField()[patchi];
-
-            const scalarField& muw = mu().boundaryField()[patchi];
-            const scalarField& mutw = mut_.boundaryField()[patchi];
-
-            scalarField magFaceGradU =
-                mag(U_.boundaryField()[patchi].snGrad());
-
-            forAll(curPatch, facei)
-            {
-                label faceCelli = curPatch.faceCells()[facei];
-
-                scalar yPlus =
-                    Cmu25*y_[faceCelli]
-                    *sqrt(k_[faceCelli])
-                    /(muw[facei]/rhow[facei]);
-
-                // For corner cells (with two boundary or more faces),
-                // omega and G in the near-wall cell are calculated
-                // as an average
-
-                cellBoundaryFaceCount[faceCelli]++;
-
-                omega_[faceCelli] +=
-                    sqrt(k_[faceCelli])
-                   /(Cmu25*kappa_.value()*y_[faceCelli]);
-
-                if (yPlus > yPlusLam_)
-                {
-                    G[faceCelli] +=
-                        (mutw[facei] + muw[facei])
-                       *magFaceGradU[facei]
-                       *Cmu25*sqrt(k_[faceCelli])
-                       /(kappa_.value()*y_[faceCelli]);
-                }
-            }
-        }
-    }
-
-
-    // Perform the averaging
-
-    forAll(patches, patchi)
-    {
-        const fvPatch& curPatch = patches[patchi];
-
-        if (isType<wallFvPatch>(curPatch))
-        {
-            forAll(curPatch, facei)
-            {
-                label faceCelli = curPatch.faceCells()[facei];
-
-                omega_[faceCelli] /= cellBoundaryFaceCount[faceCelli];
-                G[faceCelli] /= cellBoundaryFaceCount[faceCelli];
-            }
-        }
-    }
-}
-
-
-// ************************************************************************* //

src/turbulenceModels/RAS/compressible/kOmegaSST/kOmegaWallViscosityI.H

@@ -1,73 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Global
-    kOmegaWallViscosity
-
-Description
-    Calculate wall viscosity from wall-functions.
-
-\*---------------------------------------------------------------------------*/
-
-{
-    scalar Cmu25 = pow(Cmu_.value(), 0.25);
-
-    const fvPatchList& patches = mesh_.boundary();
-
-    forAll(patches, patchi)
-    {
-        const fvPatch& curPatch = patches[patchi];
-
-        if (isType<wallFvPatch>(curPatch))
-        {
-            const scalarField& rhow = rho_.boundaryField()[patchi];
-
-            const scalarField& muw = mu().boundaryField()[patchi];
-            scalarField& mutw = mut_.boundaryField()[patchi];
-
-            forAll(curPatch, facei)
-            {
-                label faceCelli = curPatch.faceCells()[facei];
-
-                scalar yPlus =
-                    Cmu25*y_[faceCelli]
-                   *sqrt(k_[faceCelli])/(muw[facei]/rhow[facei]);
-
-                if (yPlus > yPlusLam_)
-                {
-                    mutw[facei] =
-                         muw[facei]
-                        *(yPlus*kappa_.value()/log(E_.value()*yPlus) - 1);
-                }
-                else
-                {
-                    mutw[facei] = 0.0;
-                }
-            }
-        }
-    }
-}
-
-
-// ************************************************************************* //

src/turbulenceModels/RAS/compressible/realizableKE/realizableKE.C

@@ -28,6 +28,8 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "wallFvPatch.H"
 
+#include "backwardsCompatibilityWallFunctions.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -157,12 +159,11 @@ realizableKE::realizableKE
             "k",
             runTime_.timeName(),
             mesh_,
-            IOobject::MUST_READ,
+            IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh_
+        autoCreateK("k", mesh_)
     ),
-
     epsilon_
     (
         IOobject
@@ -170,12 +171,11 @@ realizableKE::realizableKE
             "epsilon",
             runTime_.timeName(),
             mesh_,
-            IOobject::MUST_READ,
+            IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh_
+        autoCreateEpsilon("epsilon", mesh_)
     ),
-
     mut_
     (
         IOobject
@@ -184,14 +184,16 @@ realizableKE::realizableKE
             runTime_.timeName(),
             mesh_,
             IOobject::NO_READ,
-            IOobject::NO_WRITE
+            IOobject::AUTO_WRITE
         ),
-        rCmu(fvc::grad(U_))*rho_*sqr(k_)/(epsilon_ + epsilonSmall_)
+        autoCreateMut("mut", mesh_)
     )
 {
     bound(k_, k0_);
     bound(epsilon_, epsilon0_);
-#   include "wallViscosityI.H"
+
+    mut_ = rCmu(fvc::grad(U_))*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    mut_.correctBoundaryConditions();
 
     printCoeffs();
 }
@@ -275,6 +277,7 @@ void realizableKE::correct()
     {
         // Re-calculate viscosity
         mut_ = rCmu(fvc::grad(U_))*rho_*sqr(k_)/epsilon_;
+        mut_.correctBoundaryConditions();
         return;
     }
 
@@ -294,9 +297,10 @@ void realizableKE::correct()
     volScalarField eta = magS*k_/epsilon_;
     volScalarField C1 = max(eta/(scalar(5) + eta), scalar(0.43));
 
-    volScalarField G = mut_*(gradU && dev(twoSymm(gradU)));
+    volScalarField G("G", mut_*(gradU && dev(twoSymm(gradU))));
 
-#   include "wallFunctionsI.H"
+    // Update espsilon and G at the wall
+    epsilon_.boundaryField().updateCoeffs();
 
     // Dissipation equation
     tmp<fvScalarMatrix> epsEqn
@@ -315,7 +319,7 @@ void realizableKE::correct()
 
     epsEqn().relax();
 
-#   include "wallDissipationI.H"
+    epsEqn().boundaryManipulate(epsilon_.boundaryField());
 
     solve(epsEqn);
     bound(epsilon_, epsilon0_);
@@ -339,9 +343,7 @@ void realizableKE::correct()
 
     // Re-calculate viscosity
     mut_ = rCmu(gradU, S2, magS)*rho_*sqr(k_)/epsilon_;
-
+    mut_.correctBoundaryConditions();
-#   include "wallViscosityI.H"
-
 }
 
 

src/turbulenceModels/RAS/incompressible/LRR/LRR.C

@@ -28,6 +28,8 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "wallFvPatch.H"
 
+#include "backwardsCompatibilityWallFunctions.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -143,12 +145,11 @@ LRR::LRR
             "R",
             runTime_.timeName(),
             mesh_,
-            IOobject::MUST_READ,
+            IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh_
+        autoCreateR("R", mesh_)
     ),
-
     k_
     (
         IOobject
@@ -156,12 +157,11 @@ LRR::LRR
             "k",
             runTime_.timeName(),
             mesh_,
-            IOobject::MUST_READ,
+            IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh_
+        autoCreateK("k", mesh_)
     ),
-
     epsilon_
     (
         IOobject
@@ -169,15 +169,26 @@ LRR::LRR
             "epsilon",
             runTime_.timeName(),
             mesh_,
-            IOobject::MUST_READ,
+            IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh_
+        autoCreateEpsilon("epsilon", mesh_)
     ),
-
+    nut_
-    nut_(Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_))
+    (
+        IOobject
+        (
+            "nut",
+            runTime_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        autoCreateNut("nut", mesh_)
+    )
 {
-#   include "wallViscosityI.H"
+    nut_ = Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    nut_.correctBoundaryConditions();
 
     if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
     {
@@ -289,9 +300,10 @@ void LRR::correct()
     RASModel::correct();
 
     volSymmTensorField P = -twoSymm(R_ & fvc::grad(U_));
-    volScalarField G = 0.5*tr(P);
+    volScalarField G("G", 0.5*tr(P));
 
-#   include "wallFunctionsI.H"
+    // Update espsilon and G at the wall
+    epsilon_.boundaryField().updateCoeffs();
 
     // Dissipation equation
     tmp<fvScalarMatrix> epsEqn
@@ -307,7 +319,7 @@ void LRR::correct()
 
     epsEqn().relax();
 
-#   include "wallDissipationI.H"
+    epsEqn().boundaryManipulate(epsilon_.boundaryField());
 
     solve(epsEqn);
     bound(epsilon_, epsilon0_);
@@ -370,8 +382,7 @@ void LRR::correct()
 
     // Re-calculate viscosity
     nut_ = Cmu_*sqr(k_)/epsilon_;
-
+    nut_.correctBoundaryConditions();
-#   include "wallViscosityI.H"
 
 
     // Correct wall shear stresses

src/turbulenceModels/RAS/incompressible/LaunderGibsonRSTM/LaunderGibsonRSTM.C

@@ -28,6 +28,8 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "wallFvPatch.H"
 
+#include "backwardsCompatibilityWallFunctions.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -172,12 +174,11 @@ LaunderGibsonRSTM::LaunderGibsonRSTM
             "R",
             runTime_.timeName(),
             mesh_,
-            IOobject::MUST_READ,
+            IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh_
+        autoCreateR("R", mesh_)
     ),
-
     k_
     (
         IOobject
@@ -185,12 +186,11 @@ LaunderGibsonRSTM::LaunderGibsonRSTM
             "k",
             runTime_.timeName(),
             mesh_,
-            IOobject::MUST_READ,
+            IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh_
+        autoCreateK("k", mesh_)
     ),
-
     epsilon_
     (
         IOobject
@@ -198,15 +198,26 @@ LaunderGibsonRSTM::LaunderGibsonRSTM
             "epsilon",
             runTime_.timeName(),
             mesh_,
-            IOobject::MUST_READ,
+            IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh_
+        autoCreateEpsilon("epsilon", mesh_)
     ),
-
+    nut_
-    nut_(Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_))
+    (
+        IOobject
+        (
+            "nut",
+            runTime_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        autoCreateNut("nut", mesh_)
+    )
 {
-#   include "wallViscosityI.H"
+    nut_ = Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    nut_.correctBoundaryConditions();
 
     if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
     {
@@ -321,9 +332,10 @@ void LaunderGibsonRSTM::correct()
     }
 
     volSymmTensorField P = -twoSymm(R_ & fvc::grad(U_));
-    volScalarField G = 0.5*tr(P);
+    volScalarField G("G", 0.5*tr(P));
 
-#   include "wallFunctionsI.H"
+    // Update espsilon and G at the wall
+    epsilon_.boundaryField().updateCoeffs();
 
     // Dissipation equation
     tmp<fvScalarMatrix> epsEqn
@@ -339,7 +351,7 @@ void LaunderGibsonRSTM::correct()
 
     epsEqn().relax();
 
-#   include "wallDissipationI.H"
+    epsEqn().boundaryManipulate(epsilon_.boundaryField());
 
     solve(epsEqn);
     bound(epsilon_, epsilon0_);
@@ -411,9 +423,7 @@ void LaunderGibsonRSTM::correct()
 
     // Re-calculate turbulent viscosity
     nut_ = Cmu_*sqr(k_)/epsilon_;
-
+    nut_.correctBoundaryConditions();
-
-#   include "wallViscosityI.H"
 
 
     // Correct wall shear stresses

src/turbulenceModels/RAS/incompressible/LienCubicKE/LienCubicKE.C

@@ -28,6 +28,8 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "wallFvPatch.H"
 
+#include "backwardsCompatibilityWallFunctions.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -152,12 +154,11 @@ LienCubicKE::LienCubicKE
             "k",
             runTime_.timeName(),
             mesh_,
-            IOobject::MUST_READ,
+            IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh_
+        autoCreateK("k", mesh_)
     ),
-
     epsilon_
     (
         IOobject
@@ -165,10 +166,10 @@ LienCubicKE::LienCubicKE
             "epsilon",
             runTime_.timeName(),
             mesh_,
-            IOobject::MUST_READ,
+            IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh_
+        autoCreateEpsilon("epsilon", mesh_)
     ),
 
     gradU_(fvc::grad(U)),
@@ -179,13 +180,25 @@ LienCubicKE::LienCubicKE
 
     C5viscosity_
     (
-        -2.0*pow(Cmu_, 3.0)*pow(k_, 4.0)/pow(epsilon_, 3.0)*
+      - 2.0*pow3(Cmu_)*pow4(k_)/pow3(epsilon_)
-        (magSqr(gradU_ + gradU_.T()) - magSqr(gradU_ - gradU_.T()))
+       *(
+            magSqr(gradU_ + gradU_.T())
+          - magSqr(gradU_ - gradU_.T())
+        )
     ),
 
-    // C5 term, implicit
+    nut_
-    nut_(Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_) + C5viscosity_),
+    (
-    // turbulent viscosity, with implicit part of C5
+        IOobject
+        (
+            "nut",
+            runTime_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        autoCreateNut("nut", mesh_)
+    ),
 
     nonlinearStress_
     (
@@ -215,7 +228,8 @@ LienCubicKE::LienCubicKE
         )
     )
 {
-#   include "wallNonlinearViscosityI.H"
+    nut_ = Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_) + C5viscosity_;
+    nut_.correctBoundaryConditions();
 
     printCoeffs();
 }
@@ -315,9 +329,14 @@ void LienCubicKE::correct()
     // generation term
     volScalarField S2 = symm(gradU_) && gradU_;
 
-    volScalarField G = Cmu_*sqr(k_)/epsilon_*S2 - (nonlinearStress_ && gradU_);
+    volScalarField G
+    (
+        "G",
+        Cmu_*sqr(k_)/epsilon_*S2 - (nonlinearStress_ && gradU_)
+    );
 
-#   include "nonLinearWallFunctionsI.H"
+    // Update espsilon and G at the wall
+    epsilon_.boundaryField().updateCoeffs();
 
     // Dissipation equation
     tmp<fvScalarMatrix> epsEqn
@@ -332,7 +351,7 @@ void LienCubicKE::correct()
 
     epsEqn().relax();
 
-#   include "wallDissipationI.H"
+    epsEqn().boundaryManipulate(epsilon_.boundaryField());
 
     solve(epsEqn);
     bound(epsilon_, epsilon0_);
@@ -367,8 +386,7 @@ void LienCubicKE::correct()
        *(magSqr(gradU_ + gradU_.T()) - magSqr(gradU_ - gradU_.T()));
 
     nut_ = Cmu_*sqr(k_)/epsilon_ + C5viscosity_;
-
+    nut_.correctBoundaryConditions();
-#   include "wallNonlinearViscosityI.H"
 
     nonlinearStress_ = symm
     (