Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

applications/solvers/combustion/fireFoam/pEqn.H

@@ -4,6 +4,9 @@ volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
+phi.boundaryField() =
+    fvc::interpolate(rho.boundaryField())
+   *(fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField());
 
 surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 

applications/solvers/multiphase/cavitatingFoam/UEqn.H

@@ -1,17 +1,8 @@
-    surfaceScalarField muEff
-    (
-        "muEff",
-        twoPhaseProperties.muf()
-      + fvc::interpolate(rho*turbulence->nut())
-    );
-
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
       + fvm::div(phi, U)
-      - fvm::laplacian(muEff, U)
+      + turbulence->divDevRhoReff(rho, U)
-    //- (fvc::grad(U) & fvc::grad(muf))
-      - fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
     );
 
     UEqn.relax();

applications/solvers/multiphase/compressibleInterFoam/UEqn.H

@@ -1,17 +1,8 @@
-    surfaceScalarField muEff
-    (
-        "muEff",
-        twoPhaseProperties.muf()
-      + fvc::interpolate(rho*turbulence->nut())
-    );
-
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
       + fvm::div(rhoPhi, U)
-      - fvm::laplacian(muEff, U)
+      + turbulence->divDevRhoReff(rho, U)
-      - (fvc::grad(U) & fvc::grad(muEff))
-    //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
     );
 
     UEqn.relax();

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H

@@ -22,27 +22,29 @@
     volVectorField HbyA2("HbyA2", U2);
     HbyA2 = rAU2*U2Eqn.H();
 
+    surfaceScalarField ppDrag("ppDrag", 0.0*phi1);
+
+    if (g0.value() > 0.0)
+    {
+        ppDrag -= ppMagf*fvc::snGrad(alpha1)*mesh.magSf();
+    }
+
+    if (kineticTheory.on())
+    {
+        ppDrag -=
+            fvc::interpolate(1.0/rho1)*rAlphaAU1f
+           *fvc::snGrad(kineticTheory.pa())*mesh.magSf();
+    }
+
     surfaceScalarField phiHbyA1
     (
         "phiHbyA1",
         (fvc::interpolate(HbyA1) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU1, alpha1, U1, phi1)
       + fvc::interpolate((1.0/rho1)*rAU1*dragCoeff)*phi2
+      + ppDrag
       + rAlphaAU1f*(g & mesh.Sf())
     );
-
-    if (g0.value() > 0.0)
-    {
-        phiHbyA1 -= ppMagf*fvc::snGrad(alpha1)*mesh.magSf();
-    }
-
-    if (kineticTheory.on())
-    {
-        phiHbyA1 -=
-            fvc::interpolate((1.0/rho1)*rAU1)
-           *fvc::snGrad(kineticTheory.pa())*mesh.magSf();
-    }
-
     mrfZones.relativeFlux(phiHbyA1);
 
     surfaceScalarField phiHbyA2
@@ -135,7 +137,8 @@
             U1 = HbyA1
                + fvc::reconstruct
                  (
-                     rAlphaAU1f
+                     ppDrag
+                   + rAlphaAU1f
                     *(
                          (g & mesh.Sf())
                        + mSfGradp/fvc::interpolate(rho1)

applications/solvers/multiphase/interFoam/UEqn.H

@@ -1,17 +1,8 @@
-    surfaceScalarField muEff
-    (
-        "muEff",
-        twoPhaseProperties.muf()
-      + fvc::interpolate(rho*turbulence->nut())
-    );
-
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
       + fvm::div(rhoPhi, U)
-      - fvm::laplacian(muEff, U)
+      + turbulence->divDevRhoReff(rho, U)
-      - (fvc::grad(U) & fvc::grad(muEff))
-    //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
     );
 
     UEqn.relax();

applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H

@@ -4,6 +4,8 @@
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
+    phi.boundaryField() =
+        fvc::interpolate(U.boundaryField()) & mesh.Sf().boundaryField();
 
     surfaceScalarField phiHbyA
     (

applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H

@@ -1,18 +1,9 @@
-    surfaceScalarField muEff
-    (
-        "muEff",
-        twoPhaseProperties->muf()
-      + fvc::interpolate(rho*turbulence->nut())
-    );
-
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
       + fvm::div(rhoPhi, U)
       - fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), U)
-      - fvm::laplacian(muEff, U)
+      + turbulence->divDevRhoReff(rho, U)
-      - (fvc::grad(U) & fvc::grad(muEff))
-    //- fvc::div(muEff*(fvc::interpolate(dev2(fvc::grad(U))) & mesh.Sf()))
     );
 
     UEqn.relax();

applications/solvers/multiphase/multiphaseInterFoam/UEqn.H

@@ -1,17 +1,8 @@
-    surfaceScalarField muEff
-    (
-        "muEff",
-        mixture.muf()
-      + fvc::interpolate(rho*turbulence->nut())
-    );
-
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
       + fvm::div(mixture.rhoPhi(), U)
-      - fvm::laplacian(muEff, U)
+      + turbulence->divDevRhoReff(rho, U)
-      - (fvc::grad(U) & fvc::grad(muEff))
-    //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
     );
 
     UEqn.relax();

applications/solvers/multiphase/settlingFoam/UEqn.H

@@ -11,7 +11,6 @@
         )
       - fvm::laplacian(muEff, U, "laplacian(muEff,U)")
       - (fvc::grad(U) & fvc::grad(muEff))
-    //- fvc::div(muEff*dev2(T(fvc::grad(U))))
     );
 
     UEqn.relax();

applications/solvers/multiphase/twoLiquidMixingFoam/UEqn.H

@@ -1,17 +1,8 @@
-    surfaceScalarField muEff
-    (
-        "muEff",
-        twoPhaseProperties.muf()
-      + fvc::interpolate(rho*turbulence->nut())
-    );
-
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
       + fvm::div(rhoPhi, U)
-      - fvm::laplacian(muEff, U)
+      + turbulence->divDevRhoReff(rho, U)
-      - (fvc::grad(U) & fvc::grad(muEff))
-    //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
     );
 
     UEqn.relax();

applications/utilities/mesh/manipulation/createPatch/createPatchDict

@@ -20,6 +20,7 @@ FoamFile
 // - always: order faces on coupled patches such that they are opposite. This
 //   is done for all coupled faces, not just for any patches created.
 // - optional: synchronise points on coupled patches.
+// - always: remove zero-sized (non-coupled) patches (that were not added)
 
 // 1. Create cyclic:
 // - specify where the faces should come from

applications/utilities/mesh/manipulation/polyDualMesh/polyDualMeshApp.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,7 +56,6 @@ Usage
 
 #include "argList.H"
 #include "Time.H"
-#include "timeSelector.H"
 #include "fvMesh.H"
 #include "unitConversion.H"
 #include "polyTopoChange.H"
@@ -354,7 +353,6 @@ int main(int argc, char *argv[])
 {
 #   include "addOverwriteOption.H"
     argList::noParallel();
-    timeSelector::addOptions(true, false);
 
     argList::validArgs.append("featureAngle [0-180]");
     argList::addBoolOption
@@ -376,9 +374,6 @@ int main(int argc, char *argv[])
 
 #   include "setRootCase.H"
 #   include "createTime.H"
-
-    instantList timeDirs = timeSelector::select0(runTime, args);
-
 #   include "createMesh.H"
 
     const word oldInstance = mesh.pointsInstance();

applications/utilities/preProcessing/setFields/setFields.C

@@ -27,7 +27,6 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
-#include "timeSelector.H"
 #include "Time.H"
 #include "fvMesh.H"
 #include "topoSetSource.H"
@@ -360,14 +359,8 @@ public:
 
 int main(int argc, char *argv[])
 {
-    timeSelector::addOptions();
-
 #   include "setRootCase.H"
 #   include "createTime.H"
-
-    // Get times list
-    instantList timeDirs = timeSelector::select0(runTime, args);
-
 #   include "createMesh.H"
 
     Info<< "Reading setFieldsDict\n" << endl;

applications/utilities/surface/surfaceMeshConvert/surfaceMeshConvert.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -58,7 +58,6 @@ Note
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
-#include "timeSelector.H"
 #include "Time.H"
 
 #include "MeshedSurfaces.H"

applications/utilities/surface/surfaceMeshConvertTesting/surfaceMeshConvertTesting.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -52,7 +52,6 @@ Note
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
-#include "timeSelector.H"
 #include "Time.H"
 #include "polyMesh.H"
 #include "triSurface.H"

applications/utilities/surface/surfaceMeshExport/surfaceMeshExport.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -59,7 +59,6 @@ Note
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
-#include "timeSelector.H"
 #include "Time.H"
 
 #include "MeshedSurfaces.H"

applications/utilities/surface/surfaceMeshImport/surfaceMeshImport.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -59,7 +59,6 @@ Note
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
-#include "timeSelector.H"
 #include "Time.H"
 
 #include "MeshedSurfaces.H"

applications/utilities/surface/surfaceMeshInfo/surfaceMeshInfo.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,7 +56,6 @@ Note
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
-#include "timeSelector.H"
 #include "Time.H"
 
 #include "UnsortedMeshedSurfaces.H"

src/OpenFOAM/db/IOobjects/IOdictionary/IOdictionaryIO.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -102,7 +102,7 @@ void Foam::IOdictionary::readFile(const bool masterOnly)
                 myComm.above(),
                 0,
                 Pstream::msgType(),
-                IOstream::ASCII
+                IOstream::BINARY
             );
             IOdictionary::readData(fromAbove);
         }
@@ -121,7 +121,7 @@ void Foam::IOdictionary::readFile(const bool masterOnly)
                 myComm.below()[belowI],
                 0,
                 Pstream::msgType(),
-                IOstream::ASCII
+                IOstream::BINARY
             );
             IOdictionary::writeData(toBelow);
         }

src/OpenFOAM/db/Time/timeSelector.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -147,7 +147,7 @@ void Foam::timeSelector::addOptions
     (
         "time",
         "ranges",
-        "comma-separated time ranges - eg, ':10,20,40-70,1000:'"
+        "comma-separated time ranges - eg, ':10,20,40:70,1000:'"
     );
 }
 

src/OpenFOAM/db/dictionary/dictionaryIO.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -80,6 +80,12 @@ Foam::autoPtr<Foam::dictionary> Foam::dictionary::New(Istream& is)
 
 bool Foam::dictionary::read(Istream& is, const bool keepHeader)
 {
+    // Check for empty dictionary
+    if (is.eof())
+    {
+        return true;
+    }
+
     if (!is.good())
     {
         FatalIOErrorIn("dictionary::read(Istream&, bool)", is)

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H

@@ -323,8 +323,6 @@ class globalMeshData
             // its own master. Maybe store as well?
 
             void calcGlobalCoPointSlaves() const;
-            const labelListList& globalCoPointSlaves() const;
-            const mapDistribute& globalCoPointSlavesMap() const;
 
 
         //- Disallow default bitwise copy construct
@@ -547,6 +545,11 @@ public:
                 //- Is my edge same orientation as master edge
                 const PackedBoolList& globalEdgeOrientation() const;
 
+            // Collocated point to collocated point
+
+                const labelListList& globalCoPointSlaves() const;
+                const mapDistribute& globalCoPointSlavesMap() const;
+
             // Coupled point to boundary faces. These are uncoupled boundary
             // faces only but include empty patches.
 

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalPoints.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -339,6 +339,31 @@ bool Foam::globalPoints::storeInitialInfo
 }
 
 
+void Foam::globalPoints::printProcPoint
+(
+    const labelList& patchToMeshPoint,
+    const labelPair& pointInfo
+) const
+{
+    label procI = globalIndexAndTransform::processor(pointInfo);
+    label index = globalIndexAndTransform::index(pointInfo);
+    label trafoI = globalIndexAndTransform::transformIndex(pointInfo);
+
+    Pout<< "    proc:" << procI;
+    Pout<< " localpoint:";
+    Pout<< index;
+    Pout<< " through transform:"
+        << trafoI << " bits:"
+        << globalTransforms_.decodeTransformIndex(trafoI);
+
+    if (procI == Pstream::myProcNo())
+    {
+        label meshPointI = localToMeshPoint(patchToMeshPoint, index);
+        Pout<< " at:" <<  mesh_.points()[meshPointI];
+    }
+}
+
+
 void Foam::globalPoints::printProcPoints
 (
     const labelList& patchToMeshPoint,
@@ -347,23 +372,7 @@ void Foam::globalPoints::printProcPoints
 {
     forAll(pointInfo, i)
     {
-        label procI = globalIndexAndTransform::processor(pointInfo[i]);
+        printProcPoint(patchToMeshPoint, pointInfo[i]);
-        label index = globalIndexAndTransform::index(pointInfo[i]);
-        label trafoI = globalIndexAndTransform::transformIndex(pointInfo[i]);
-
-        Pout<< "    proc:" << procI;
-        Pout<< " localpoint:";
-        Pout<< index;
-        Pout<< " through transform:"
-            << trafoI << " bits:"
-            << globalTransforms_.decodeTransformIndex(trafoI);
-
-        if (procI == Pstream::myProcNo())
-        {
-            label meshPointI = localToMeshPoint(patchToMeshPoint, index);
-            Pout<< " at:" <<  mesh_.points()[meshPointI];
-        }
-
         Pout<< endl;
     }
 }

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalPoints.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -206,6 +206,12 @@ class globalPoints
         );
 
         //- Debug printing
+        void printProcPoint
+        (
+            const labelList& patchToMeshPoint,
+            const labelPair& pointInfo
+        ) const;
+
         void printProcPoints
         (
             const labelList& patchToMeshPoint,

src/Pstream/mpi/UIPread.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -107,12 +107,7 @@ Foam::UIPstream::UIPstream
 
         if (!messageSize_)
         {
-            FatalErrorIn
+            setEof();
-            (
-                "UIPstream::UIPstream(const commsTypes, const int, "
-                "DynamicList<char>&, streamFormat, versionNumber)"
-            )   << "read failed"
-                << Foam::abort(FatalError);
         }
     }
 }
@@ -199,11 +194,7 @@ Foam::UIPstream::UIPstream(const int fromProcNo, PstreamBuffers& buffers)
 
         if (!messageSize_)
         {
-            FatalErrorIn
+            setEof();
-            (
-                "UIPstream::UIPstream(const int, PstreamBuffers&)"
-            )   << "read failed"
-                << Foam::abort(FatalError);
         }
     }
 }

src/finiteVolume/fields/fvPatchFields/basic/calculated/calculatedFvPatchField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -98,20 +98,19 @@ calculatedFvPatchField<Type>::calculatedFvPatchField
 
 
 template<class Type>
-template<class Type2>
 tmp<fvPatchField<Type> > fvPatchField<Type>::NewCalculatedType
 (
-    const fvPatchField<Type2>& pf
+    const fvPatch& p
 )
 {
     typename patchConstructorTable::iterator patchTypeCstrIter =
-        patchConstructorTablePtr_->find(pf.patch().type());
+        patchConstructorTablePtr_->find(p.type());
 
     if (patchTypeCstrIter != patchConstructorTablePtr_->end())
     {
         return patchTypeCstrIter()
         (
-            pf.patch(),
+            p,
             DimensionedField<Type, volMesh>::null()
         );
     }
@@ -121,7 +120,7 @@ tmp<fvPatchField<Type> > fvPatchField<Type>::NewCalculatedType
         (
             new calculatedFvPatchField<Type>
             (
-                pf.patch(),
+                p,
                 DimensionedField<Type, volMesh>::null()
             )
         );
@@ -129,6 +128,17 @@ tmp<fvPatchField<Type> > fvPatchField<Type>::NewCalculatedType
 }
 
 
+template<class Type>
+template<class Type2>
+tmp<fvPatchField<Type> > fvPatchField<Type>::NewCalculatedType
+(
+    const fvPatchField<Type2>& pf
+)
+{
+    return NewCalculatedType(pf.patch());
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>

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

@@ -257,6 +257,13 @@ public:
             const dictionary&
         );
 
+        //- Return a pointer to a new calculatedFvPatchField created on
+        //  freestore without setting patchField values
+        static tmp<fvPatchField<Type> > NewCalculatedType
+        (
+            const fvPatch&
+        );
+
         //- Return a pointer to a new calculatedFvPatchField created on
         //  freestore without setting patchField values
         template<class Type2>

src/finiteVolume/fields/fvsPatchFields/basic/calculated/calculatedFvsPatchField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -97,21 +97,20 @@ calculatedFvsPatchField<Type>::calculatedFvsPatchField
 
 
 template<class Type>
-template<class Type2>
 tmp<fvsPatchField<Type> > fvsPatchField<Type>::NewCalculatedType
 (
-    const fvsPatchField<Type2>& pf
+    const fvPatch& p
 )
 {
     typename patchConstructorTable::iterator patchTypeCstrIter =
-        patchConstructorTablePtr_->find(pf.patch().type());
+        patchConstructorTablePtr_->find(p.type());
 
     if (patchTypeCstrIter != patchConstructorTablePtr_->end())
     {
         return patchTypeCstrIter()
         (
-            pf.patch(),
+            p,
-            Field<Type>::null()
+            DimensionedField<Type, surfaceMesh>::null()
         );
     }
     else
@@ -120,14 +119,25 @@ tmp<fvsPatchField<Type> > fvsPatchField<Type>::NewCalculatedType
         (
             new calculatedFvsPatchField<Type>
             (
-                pf.patch(),
+                p,
-                Field<Type>::null()
+                DimensionedField<Type, surfaceMesh>::null()
             )
         );
     }
 }
 
 
+template<class Type>
+template<class Type2>
+tmp<fvsPatchField<Type> > fvsPatchField<Type>::NewCalculatedType
+(
+    const fvsPatchField<Type2>& pf
+)
+{
+    return NewCalculatedType(pf.patch());
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/finiteVolume/fields/fvsPatchFields/fvsPatchField/fvsPatchField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -247,6 +247,13 @@ public:
             const dictionary&
         );
 
+        //- Return a pointer to a new calculatedFvsPatchField created on
+        //  freestore without setting patchField values
+        static tmp<fvsPatchField<Type> > NewCalculatedType
+        (
+            const fvPatch&
+        );
+
         //- Return a pointer to a new calculatedFvsPatchField created on
         //  freestore without setting patchField values
         template<class Type2>

src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -302,6 +302,43 @@ interpolate
 }
 
 
+template<class Type>
+tmp<FieldField<fvsPatchField, Type> > interpolate
+(
+    const FieldField<fvPatchField, Type>& fvpff
+)
+{
+    FieldField<fvsPatchField, Type>* fvspffPtr
+    (
+        new FieldField<fvsPatchField, Type>(fvpff.size())
+    );
+
+    forAll(*fvspffPtr, patchi)
+    {
+        fvspffPtr->set
+        (
+            patchi,
+            fvsPatchField<Type>::NewCalculatedType(fvpff[patchi].patch()).ptr()
+        );
+        (*fvspffPtr)[patchi] = fvpff[patchi];
+    }
+
+    return tmp<FieldField<fvsPatchField, Type> >(fvspffPtr);
+}
+
+
+template<class Type>
+tmp<FieldField<fvsPatchField, Type> > interpolate
+(
+    const tmp<FieldField<fvPatchField, Type> >& tfvpff
+)
+{
+    tmp<FieldField<fvPatchField, Type> > tfvspff = interpolate(tfvpff());
+    tfvpff.clear();
+    return tfvspff;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace fvc

src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolate.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -163,12 +163,26 @@ namespace fvc
         const tmp<GeometricField<Type, fvPatchField, volMesh> >& tvf
     );
 
-    //- Interpolate tmp field onto faces using 'interpolate(\<name\>)'
+    //- Interpolate field onto faces using 'interpolate(\<name\>)'
     template<class Type>
     static tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > interpolate
     (
         const GeometricField<Type, fvPatchField, volMesh>& tvf
     );
+
+    //- Interpolate boundary field onto faces (simply a type conversion)
+    template<class Type>
+    static tmp<FieldField<fvsPatchField, Type> > interpolate
+    (
+        const FieldField<fvPatchField, Type>& fvpff
+    );
+
+    //- Interpolate boundary field onto faces (simply a type conversion)
+    template<class Type>
+    static tmp<FieldField<fvsPatchField, Type> > interpolate
+    (
+        const tmp<FieldField<fvPatchField, Type> >& tfvpff
+    );
 }
 
 

src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolate.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -51,8 +51,8 @@ void volPointInterpolation::syncUntransformedData
     const indirectPrimitivePatch& cpp = gmd.coupledPatch();
     const labelList& meshPoints = cpp.meshPoints();
 
-    const mapDistribute& slavesMap = gmd.globalPointSlavesMap();
+    const mapDistribute& slavesMap = gmd.globalCoPointSlavesMap();
-    const labelListList& slaves = gmd.globalPointSlaves();
+    const labelListList& slaves = gmd.globalCoPointSlaves();
 
     List<Type> elems(slavesMap.constructSize());
     forAll(meshPoints, i)
@@ -105,8 +105,8 @@ void volPointInterpolation::pushUntransformedData
     const indirectPrimitivePatch& cpp = gmd.coupledPatch();
     const labelList& meshPoints = cpp.meshPoints();
 
-    const mapDistribute& slavesMap = gmd.globalPointSlavesMap();
+    const mapDistribute& slavesMap = gmd.globalCoPointSlavesMap();
-    const labelListList& slaves = gmd.globalPointSlaves();
+    const labelListList& slaves = gmd.globalCoPointSlaves();
 
     List<Type> elems(slavesMap.constructSize());
     forAll(meshPoints, i)
@@ -155,14 +155,14 @@ void volPointInterpolation::addSeparated
             refCast<coupledPointPatchField<Type> >
                 (pf.boundaryField()[patchI]).initSwapAddSeparated
                 (
-                    Pstream::blocking,  //Pstream::nonBlocking,
+                    Pstream::nonBlocking,
                     pf.internalField()
                 );
         }
     }
 
     // Block for any outstanding requests
-    //Pstream::waitRequests();
+    Pstream::waitRequests();
 
     forAll(pf.boundaryField(), patchI)
     {
@@ -171,7 +171,7 @@ void volPointInterpolation::addSeparated
             refCast<coupledPointPatchField<Type> >
                 (pf.boundaryField()[patchI]).swapAddSeparated
                 (
-                    Pstream::blocking,  //Pstream::nonBlocking,
+                    Pstream::nonBlocking,
                     pf.internalField()
                 );
         }
@@ -306,7 +306,6 @@ void volPointInterpolation::interpolateBoundaryField
     }
 
     // Sum collocated contributions
-    //mesh().globalData().syncPointData(pfi, plusEqOp<Type>());
     syncUntransformedData(pfi, plusEqOp<Type>());
 
     // And add separated contributions
@@ -314,9 +313,7 @@ void volPointInterpolation::interpolateBoundaryField
 
     // Push master data to slaves. It is possible (not sure how often) for
     // a coupled point to have its master on a different patch so
-    // to make sure just push master data to slaves. Reuse the syncPointData
+    // to make sure just push master data to slaves.
-    // structure.
-    //mesh().globalData().syncPointData(pfi, nopEqOp<Type>());
     pushUntransformedData(pfi);
 
 

src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolation.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -113,10 +113,8 @@ void volPointInterpolation::calcBoundaryAddressing()
     {
         boolList oldData(isPatchPoint_);
 
-        //mesh().globalData().syncPointData(isPatchPoint_, orEqOp<bool>());
         syncUntransformedData(isPatchPoint_, orEqOp<bool>());
 
-
         forAll(isPatchPoint_, pointI)
         {
             if (isPatchPoint_[pointI] != oldData[pointI])
@@ -272,7 +270,7 @@ void volPointInterpolation::makeWeights()
     makeInternalWeights(sumWeights);
 
 
-    // Create boundary weights; add to sumWeights
+    // Create boundary weights; override sumWeights
     makeBoundaryWeights(sumWeights);
 
 
@@ -292,7 +290,6 @@ void volPointInterpolation::makeWeights()
 
 
     // Sum collocated contributions
-    //mesh().globalData().syncPointData(sumWeights, plusEqOp<scalar>());
     syncUntransformedData(sumWeights, plusEqOp<scalar>());
 
     // And add separated contributions
@@ -302,7 +299,6 @@ void volPointInterpolation::makeWeights()
     // a coupled point to have its master on a different patch so
     // to make sure just push master data to slaves. Reuse the syncPointData
     // structure.
-    //mesh().globalData().syncPointData(sumWeights, nopEqOp<scalar>());
     pushUntransformedData(sumWeights);
 
 

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C

@@ -713,7 +713,7 @@ Foam::label Foam::AMIInterpolation<SourcePatch, TargetPatch>::findTargetFace
 
     const pointField& srcPts = srcPatch.points();
     const face& srcFace = srcPatch[srcFaceI];
-    const point& srcPt = srcFace.centre(srcPts);
+    const point srcPt = srcFace.centre(srcPts);
     const scalar srcFaceArea = srcMagSf_[srcFaceI];
 
 //    pointIndexHit sample = treePtr_->findNearest(srcPt, sqr(0.1*bb.mag()));
@@ -782,6 +782,62 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::appendNbrFaces
 }
 
 
+template<class SourcePatch, class TargetPatch>
+bool Foam::AMIInterpolation<SourcePatch, TargetPatch>::processSourceFace
+(
+    const SourcePatch& srcPatch,
+    const TargetPatch& tgtPatch,
+    const label srcFaceI,
+    const label tgtStartFaceI,
+
+    // list of tgt face neighbour faces
+    DynamicList<label>& nbrFaces,
+    // list of faces currently visited for srcFaceI to avoid multiple hits
+    DynamicList<label>& visitedFaces,
+
+    // temporary storage for addressing and weights
+    List<DynamicList<label> >& srcAddr,
+    List<DynamicList<scalar> >& srcWght,
+    List<DynamicList<label> >& tgtAddr,
+    List<DynamicList<scalar> >& tgtWght
+)
+{
+    nbrFaces.clear();
+    visitedFaces.clear();
+
+    // append initial target face and neighbours
+    nbrFaces.append(tgtStartFaceI);
+    appendNbrFaces(tgtStartFaceI, tgtPatch, visitedFaces, nbrFaces);
+
+    bool faceProcessed = false;
+
+    do
+    {
+        // process new target face
+        label tgtFaceI = nbrFaces.remove();
+        visitedFaces.append(tgtFaceI);
+        scalar area = interArea(srcFaceI, tgtFaceI, srcPatch, tgtPatch);
+
+        // store when intersection area > 0
+        if (area > 0)
+        {
+            srcAddr[srcFaceI].append(tgtFaceI);
+            srcWght[srcFaceI].append(area);
+
+            tgtAddr[tgtFaceI].append(srcFaceI);
+            tgtWght[tgtFaceI].append(area);
+
+            appendNbrFaces(tgtFaceI, tgtPatch, visitedFaces, nbrFaces);
+
+            faceProcessed = true;
+        }
+
+    } while (nbrFaces.size() > 0);
+
+    return faceProcessed;
+}
+
+
 template<class SourcePatch, class TargetPatch>
 void Foam::AMIInterpolation<SourcePatch, TargetPatch>::setNextFaces
 (
@@ -811,7 +867,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::setNextFaces
                 label faceT = visitedFaces[j];
                 scalar area = interArea(faceS, faceT, srcPatch0, tgtPatch0);
 
-                if (area > 0)
+                // Check that faces have enough overlap for robust walking
+                if (area/srcMagSf_[srcFaceI] > faceAreaIntersect::tolerance())
                 {
                     // TODO - throwing area away - re-use in next iteration?
 
@@ -864,7 +921,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::setNextFaces
                 << "target face" << endl;
         }
 
-//        foundNextSeed = false;
+        foundNextSeed = false;
         for (label faceI = startSeedI; faceI < mapFlag.size(); faceI++)
         {
             if (mapFlag[faceI])
@@ -887,7 +944,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::setNextFaces
 
         FatalErrorIn
         (
-            "void Foam::cyclicAMIPolyPatch::setNextFaces"
+            "void Foam::AMIInterpolation<SourcePatch, TargetPatch>::"
+            "setNextFaces"
             "("
                 "label&, "
                 "label&, "
@@ -946,6 +1004,25 @@ Foam::scalar Foam::AMIInterpolation<SourcePatch, TargetPatch>::interArea
     {
         area = inter.calc(src, tgt, n, triMode_);
     }
+    else
+    {
+        WarningIn
+        (
+            "void Foam::AMIInterpolation<SourcePatch, TargetPatch>::"
+            "interArea"
+            "("
+                "const label, "
+                "const label, "
+                "const SourcePatch&, "
+                "const TargetPatch&"
+            ") const"
+        )   << "Invalid normal for source face " << srcFaceI
+            << " points " << UIndirectList<point>(srcPoints, src)
+            << " target face " << tgtFaceI
+            << " points " << UIndirectList<point>(tgtPoints, tgt)
+            << endl;
+    }
+
 
     if ((debug > 1) && (area > 0))
     {
@@ -956,6 +1033,105 @@ Foam::scalar Foam::AMIInterpolation<SourcePatch, TargetPatch>::interArea
 }
 
 
+template<class SourcePatch, class TargetPatch>
+void Foam::AMIInterpolation<SourcePatch, TargetPatch>::
+restartUncoveredSourceFace
+(
+    const SourcePatch& srcPatch,
+    const TargetPatch& tgtPatch,
+    List<DynamicList<label> >& srcAddr,
+    List<DynamicList<scalar> >& srcWght,
+    List<DynamicList<label> >& tgtAddr,
+    List<DynamicList<scalar> >& tgtWght
+)
+{
+    // Collect all src faces with a low weight
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    labelHashSet lowWeightFaces(100);
+    forAll(srcWght, srcFaceI)
+    {
+        scalar s = sum(srcWght[srcFaceI]);
+        scalar t = s/srcMagSf_[srcFaceI];
+
+        if (t < 0.5)
+        {
+            lowWeightFaces.insert(srcFaceI);
+        }
+    }
+
+    Info<< "AMIInterpolation : restarting search on "
+        << returnReduce(lowWeightFaces.size(), sumOp<label>())
+        << " faces since sum of weights < 0.5" << endl;
+
+    if (lowWeightFaces.size() > 0)
+    {
+        // Erase all the lowWeight source faces from the target
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        DynamicList<label> okSrcFaces(10);
+        DynamicList<scalar> okSrcWeights(10);
+        forAll(tgtAddr, tgtFaceI)
+        {
+            okSrcFaces.clear();
+            okSrcWeights.clear();
+            DynamicList<label>& srcFaces = tgtAddr[tgtFaceI];
+            DynamicList<scalar>& srcWeights = tgtWght[tgtFaceI];
+            forAll(srcFaces, i)
+            {
+                if (!lowWeightFaces.found(srcFaces[i]))
+                {
+                    okSrcFaces.append(srcFaces[i]);
+                    okSrcWeights.append(srcWeights[i]);
+                }
+            }
+            if (okSrcFaces.size() < srcFaces.size())
+            {
+                srcFaces.transfer(okSrcFaces);
+                srcWeights.transfer(okSrcWeights);
+            }
+        }
+
+
+
+        // Restart search from best hit
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        // list of tgt face neighbour faces
+        DynamicList<label> nbrFaces(10);
+
+        // list of faces currently visited for srcFaceI to avoid multiple hits
+        DynamicList<label> visitedFaces(10);
+
+        forAllConstIter(labelHashSet, lowWeightFaces, iter)
+        {
+            label srcFaceI = iter.key();
+            label tgtFaceI = findTargetFace(srcFaceI, srcPatch);
+            if (tgtFaceI != -1)
+            {
+                //bool faceProcessed =
+                processSourceFace
+                (
+                    srcPatch,
+                    tgtPatch,
+                    srcFaceI,
+                    tgtFaceI,
+
+                    nbrFaces,
+                    visitedFaces,
+
+                    srcAddr,
+                    srcWght,
+                    tgtAddr,
+                    tgtWght
+                );
+                // ? Check faceProcessed to see if restarting has worked.
+            }
+        }
+    }
+}
+
+
 template<class SourcePatch, class TargetPatch>
 void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
 (
@@ -1073,37 +1249,22 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
     DynamicList<label> nonOverlapFaces;
     do
     {
-        nbrFaces.clear();
+        // Do advancing front starting from srcFaceI,tgtFaceI
-        visitedFaces.clear();
+        bool faceProcessed = processSourceFace
+        (
+            srcPatch,
+            tgtPatch,
+            srcFaceI,
+            tgtFaceI,
 
-        // append initial target face and neighbours
+            nbrFaces,
-        nbrFaces.append(tgtFaceI);
+            visitedFaces,
-        appendNbrFaces(tgtFaceI, tgtPatch, visitedFaces, nbrFaces);
 
-        bool faceProcessed = false;
+            srcAddr,
-
+            srcWght,
-        do
+            tgtAddr,
-        {
+            tgtWght
-            // process new target face
+        );
-            tgtFaceI = nbrFaces.remove();
-            visitedFaces.append(tgtFaceI);
-            scalar area = interArea(srcFaceI, tgtFaceI, srcPatch, tgtPatch);
-
-            // store when intersection area > 0
-            if (area > 0)
-            {
-                srcAddr[srcFaceI].append(tgtFaceI);
-                srcWght[srcFaceI].append(area);
-
-                tgtAddr[tgtFaceI].append(srcFaceI);
-                tgtWght[tgtFaceI].append(area);
-
-                appendNbrFaces(tgtFaceI, tgtPatch, visitedFaces, nbrFaces);
-
-                faceProcessed = true;
-            }
-
-        } while (nbrFaces.size() > 0);
 
         mapFlag[srcFaceI] = false;
 
@@ -1140,6 +1301,22 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
         srcNonOverlap_.transfer(nonOverlapFaces);
     }
 
+
+    // Check for badly covered faces
+    if (debug)
+    {
+        restartUncoveredSourceFace
+        (
+            srcPatch,
+            tgtPatch,
+            srcAddr,
+            srcWght,
+            tgtAddr,
+            tgtWght
+        );
+    }
+
+
     // transfer data to persistent storage
     forAll(srcAddr, i)
     {

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H

@@ -246,6 +246,31 @@ class AMIInterpolation
                 DynamicList<label>& faceIDs
             ) const;
 
+            bool processSourceFace
+            (
+                const SourcePatch& srcPatch,
+                const TargetPatch& tgtPatch,
+                const label srcFaceI,
+                const label tgtStartFaceI,
+
+                DynamicList<label>& nbrFaces,
+                DynamicList<label>& visitedFaces,
+                List<DynamicList<label> >& srcAddr,
+                List<DynamicList<scalar> >& srcWght,
+                List<DynamicList<label> >& tgtAddr,
+                List<DynamicList<scalar> >& tgtWght
+            );
+
+            void restartUncoveredSourceFace
+            (
+                const SourcePatch& srcPatch,
+                const TargetPatch& tgtPatch,
+                List<DynamicList<label> >& srcAddr,
+                List<DynamicList<scalar> >& srcWght,
+                List<DynamicList<label> >& tgtAddr,
+                List<DynamicList<scalar> >& tgtWght
+            );
+
             //- Set the source and target seed faces
             void setNextFaces
             (

src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,6 +25,10 @@ License
 
 #include "faceAreaIntersect.H"
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+Foam::scalar Foam::faceAreaIntersect::tol = 1e-6;
+
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
 void Foam::faceAreaIntersect::triSliceWithPlane
@@ -36,8 +40,6 @@ void Foam::faceAreaIntersect::triSliceWithPlane
     const scalar len
 )
 {
-    const scalar matchTol = 1e-6;
-
     // distance to cutting plane
     FixedList<scalar, 3> d;
 
@@ -51,7 +53,7 @@ void Foam::faceAreaIntersect::triSliceWithPlane
     {
         d[i] = ((tri[i] - p.refPoint()) & p.normal());
 
-        if (mag(d[i]) < matchTol*len)
+        if (mag(d[i]) < tol*len)
         {
             nCoPlanar++;
             copI = i;

src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -78,6 +78,11 @@ private:
         const bool reverseB_;
 
 
+    // Static data members
+
+        static scalar tol;
+
+
     // Private Member Functions
 
         //- Get triPoints from face
@@ -152,6 +157,9 @@ public:
 
     // Public Member Functions
 
+        //- Fraction of local length scale to use as intersection tolerance
+        inline static scalar& tolerance();
+
         //- Return area of intersection of faceA with faceB
         scalar calc
         (

src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,6 +23,8 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
 inline void Foam::faceAreaIntersect::setTriPoints
 (
     const point& a,
@@ -106,4 +108,12 @@ inline Foam::scalar Foam::faceAreaIntersect::triArea(const triPoints& t) const
 }
 
 
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::scalar& Foam::faceAreaIntersect::tolerance()
+{
+    return tol;
+}
+
+
 // ************************************************************************* //

src/postProcessing/foamCalcFunctions/field/magSqr/magSqr.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,7 +56,7 @@ Foam::calcTypes::magSqr::~magSqr()
 
 void Foam::calcTypes::magSqr::init()
 {
-    Foam::argList::validArgs.append("magSqr");
+    argList::validArgs.append("magSqr");
     argList::validArgs.append("fieldName");
 }
 

src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C

@@ -27,8 +27,8 @@ License
 #include "volFields.H"
 #include "regionSplit.H"
 #include "fvcVolumeIntegrate.H"
-#include "Histogram.H"
 #include "mathematicalConstants.H"
+#include "stringListOps.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -56,6 +56,281 @@ namespace Foam
 }
 
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::regionSizeDistribution::writeGraph
+(
+    const coordSet& coords,
+    const word& valueName,
+    const scalarField& values
+) const
+{
+    const fvMesh& mesh = refCast<const fvMesh>(obr_);
+
+    const wordList valNames(1, valueName);
+
+    fileName outputPath;
+    if (Pstream::parRun())
+    {
+        outputPath = mesh.time().path()/".."/name_;
+    }
+    else
+    {
+        outputPath = mesh.time().path()/name_;
+    }
+
+    if (mesh.name() != fvMesh::defaultRegion)
+    {
+        outputPath = outputPath/mesh.name();
+    }
+
+    mkDir(outputPath/mesh.time().timeName());
+    OFstream str
+    (
+        outputPath
+      / mesh.time().timeName()
+      / formatterPtr_().getFileName(coords, valNames)
+    );
+    Info<< "Writing distribution of " << valueName << " to " << str.name()
+        << endl;
+
+    List<const scalarField*> valPtrs(1);
+    valPtrs[0] = &values;
+    formatterPtr_().write(coords, valNames, valPtrs, str);
+}
+
+
+void Foam::regionSizeDistribution::writeAlphaFields
+(
+    const regionSplit& regions,
+    const Map<label>& patchRegions,
+    const Map<scalar>& regionVolume,
+    const volScalarField& alpha
+) const
+{
+    const scalar maxDropletVol = 1.0/6.0*pow(maxDiam_, 3);
+
+    // Split alpha field
+    // ~~~~~~~~~~~~~~~~~
+    // Split into
+    //  - liquidCore            : region connected to inlet patches
+    //  - per region a volume   : for all other regions
+    //  - backgroundAlpha       : remaining alpha
+
+
+    // Construct field
+    volScalarField liquidCore
+    (
+        IOobject
+        (
+            alphaName_ + "_liquidCore",
+            obr_.time().timeName(),
+            obr_,
+            IOobject::NO_READ
+        ),
+        alpha,
+        fvPatchField<scalar>::calculatedType()
+    );
+
+    volScalarField backgroundAlpha
+    (
+        IOobject
+        (
+            alphaName_ + "_background",
+            obr_.time().timeName(),
+            obr_,
+            IOobject::NO_READ
+        ),
+        alpha,
+        fvPatchField<scalar>::calculatedType()
+    );
+
+
+    // Knock out any cell not in patchRegions
+    forAll(liquidCore, cellI)
+    {
+        label regionI = regions[cellI];
+        if (patchRegions.found(regionI))
+        {
+            backgroundAlpha[cellI] = 0;
+        }
+        else
+        {
+            liquidCore[cellI] = 0;
+
+            scalar regionVol = regionVolume[regionI];
+            if (regionVol < maxDropletVol)
+            {
+                backgroundAlpha[cellI] = 0;
+            }
+        }
+    }
+    liquidCore.correctBoundaryConditions();
+    backgroundAlpha.correctBoundaryConditions();
+
+    Info<< "Volume of liquid-core = "
+        << fvc::domainIntegrate(liquidCore).value()
+        << endl;
+    Info<< "Volume of background  = "
+        << fvc::domainIntegrate(backgroundAlpha).value()
+        << endl;
+
+    Info<< "Writing liquid-core field to " << liquidCore.name() << endl;
+    liquidCore.write();
+    Info<< "Writing background field to " << backgroundAlpha.name() << endl;
+    backgroundAlpha.write();
+}
+
+
+Foam::Map<Foam::label> Foam::regionSizeDistribution::findPatchRegions
+(
+    const polyMesh& mesh,
+    const regionSplit& regions
+) const
+{
+    // Mark all regions starting at patches
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Count number of patch faces (just for initial sizing)
+    const labelHashSet patchIDs(mesh.boundaryMesh().patchSet(patchNames_));
+
+    label nPatchFaces = 0;
+    forAllConstIter(labelHashSet, patchIDs, iter)
+    {
+        nPatchFaces += mesh.boundaryMesh()[iter.key()].size();
+    }
+
+
+    Map<label> patchRegions(nPatchFaces);
+    forAllConstIter(labelHashSet, patchIDs, iter)
+    {
+        const polyPatch& pp = mesh.boundaryMesh()[iter.key()];
+
+        // Collect all regions on the patch
+        const labelList& faceCells = pp.faceCells();
+
+        forAll(faceCells, i)
+        {
+            patchRegions.insert
+            (
+                regions[faceCells[i]],
+                Pstream::myProcNo()     // dummy value
+            );
+        }
+    }
+
+
+    // Make sure all the processors have the same set of regions
+    Pstream::mapCombineGather(patchRegions, minEqOp<label>());
+    Pstream::mapCombineScatter(patchRegions);
+
+    return patchRegions;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::regionSizeDistribution::divide
+(
+    const scalarField& num,
+    const scalarField& denom
+)
+{
+    tmp<scalarField> tresult(new scalarField(num.size()));
+    scalarField& result = tresult();
+
+    forAll(denom, i)
+    {
+        if (denom[i] != 0)
+        {
+            result[i] = num[i]/denom[i];
+        }
+        else
+        {
+            result[i] = 0.0;
+        }
+    }
+    return tresult;
+}
+
+
+void Foam::regionSizeDistribution::writeGraphs
+(
+    const word& fieldName,              // name of field
+    const labelList& indices,           // index of bin for each region
+    const scalarField& sortedField,     // per region field data
+    const scalarField& binCount,        // per bin number of regions
+    const coordSet& coords              // graph data for bins
+) const
+{
+    if (Pstream::master())
+    {
+        // Calculate per-bin average
+        scalarField binSum(nBins_, 0.0);
+        forAll(sortedField, i)
+        {
+            binSum[indices[i]] += sortedField[i];
+        }
+
+        scalarField binAvg(divide(binSum, binCount));
+
+        // Per bin deviation
+        scalarField binSqrSum(nBins_, 0.0);
+        forAll(sortedField, i)
+        {
+            binSqrSum[indices[i]] += Foam::sqr(sortedField[i]);
+        }
+        scalarField binDev
+        (
+            sqrt(divide(binSqrSum, binCount) - Foam::sqr(binAvg))
+        );
+
+        // Write average
+        writeGraph(coords, fieldName + "_sum", binSum);
+        // Write average
+        writeGraph(coords, fieldName + "_avg", binAvg);
+        // Write deviation
+        writeGraph(coords, fieldName + "_dev", binDev);
+    }
+}
+
+
+void Foam::regionSizeDistribution::writeGraphs
+(
+    const word& fieldName,              // name of field
+    const scalarField& cellField,       // per cell field data
+    const regionSplit& regions,         // per cell the region(=droplet)
+    const labelList& sortedRegions,     // valid regions in sorted order
+    const scalarField& sortedNormalisation,
+
+    const labelList& indices,           // per region index of bin
+    const scalarField& binCount,        // per bin number of regions
+    const coordSet& coords              // graph data for bins
+) const
+{
+    // Sum on a per-region basis. Parallel reduced.
+    Map<scalar> regionField(regionSum(regions, cellField));
+
+    // Extract in region order
+    scalarField sortedField
+    (
+        sortedNormalisation
+      * extractData
+        (
+            sortedRegions,
+            regionField
+        )
+    );
+
+    writeGraphs
+    (
+        fieldName,      // name of field
+        indices,        // index of bin for each region
+        sortedField,    // per region field data
+        binCount,       // per bin number of regions
+        coords          // graph data for bins
+    );
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::regionSizeDistribution::regionSizeDistribution
@@ -103,8 +378,11 @@ void Foam::regionSizeDistribution::read(const dictionary& dict)
         dict.lookup("field") >> alphaName_;
         dict.lookup("patches") >> patchNames_;
         dict.lookup("threshold") >> threshold_;
-        dict.lookup("volFraction") >> volFraction_;
+        dict.lookup("maxDiameter") >> maxDiam_;
+        minDiam_ = 0.0;
+        dict.readIfPresent("minDiameter", minDiam_);
         dict.lookup("nBins") >> nBins_;
+        dict.lookup("fields") >> fields_;
 
         word format(dict.lookup("setFormat"));
         formatterPtr_ = writer<scalar>::New(format);
@@ -163,14 +441,17 @@ void Foam::regionSizeDistribution::write()
            : obr_.lookupObject<volScalarField>(alphaName_)
         );
 
-        Info<< "Volume of alpha = "
+        Info<< "Volume of alpha          = "
             << fvc::domainIntegrate(alpha).value()
             << endl;
 
         const scalar meshVol = gSum(mesh.V());
-        Info<< "Mesh volume = " << meshVol << endl;
+        const scalar maxDropletVol = 1.0/6.0*pow(maxDiam_, 3);
-        Info<< "Background region volume limit = " << volFraction_*meshVol
+        const scalar delta = (maxDiam_-minDiam_)/nBins_;
-            << endl;
+
+        Info<< "Mesh volume              = " << meshVol << endl;
+        Info<< "Maximum droplet diameter = " << maxDiam_ << endl;
+        Info<< "Maximum droplet volume   = " << maxDropletVol << endl;
 
 
         // Determine blocked faces
@@ -260,138 +541,115 @@ void Foam::regionSizeDistribution::write()
         }
 
 
+        // Determine regions connected to supplied patches
+        Map<label> patchRegions(findPatchRegions(mesh, regions));
+
+
+
         // Sum all regions
-        Map<Pair<scalar> > regionVolume(regions.nRegions()/Pstream::nProcs());
+        const scalarField alphaVol = alpha.internalField()*mesh.V();
-        forAll(alpha, cellI)
+        Map<scalar> allRegionVolume(regionSum(regions, mesh.V()));
-        {
+        Map<scalar> allRegionAlphaVolume(regionSum(regions, alphaVol));
-            scalar cellVol = mesh.V()[cellI];
+        Map<label> allRegionNumCells
-            scalar alphaVol = alpha[cellI]*cellVol;
+        (
-
+            regionSum
-            label regionI = regions[cellI];
+            (
-
+                regions,
-            Map<Pair<scalar> >::iterator fnd = regionVolume.find(regionI);
+                labelField(mesh.nCells(), 1.0)
-            if (fnd == regionVolume.end())
+            )
-            {
+        );
-                regionVolume.insert
-                (
-                    regionI,
-                    Pair<scalar>(cellVol, alphaVol)
-                );
-            }
-            else
-            {
-                fnd().first() += cellVol;
-                fnd().second() += alphaVol;
-            }
-        }
-        Pstream::mapCombineGather(regionVolume, ListPlusEqOp<scalar, 2>());
-        Pstream::mapCombineScatter(regionVolume);
-
 
         if (debug)
         {
             Info<< token::TAB << "Region"
                 << token::TAB << "Volume(mesh)"
                 << token::TAB << "Volume(" << alpha.name() << "):"
+                << token::TAB << "nCells"
                 << endl;
             scalar meshSumVol = 0.0;
             scalar alphaSumVol = 0.0;
+            label nCells = 0;
 
-            forAllConstIter(Map<Pair<scalar> >, regionVolume, iter)
+            Map<scalar>::const_iterator vIter = allRegionVolume.begin();
+            Map<scalar>::const_iterator aIter = allRegionAlphaVolume.begin();
+            Map<label>::const_iterator numIter = allRegionNumCells.begin();
+            for
+            (
+                ;
+                vIter != allRegionVolume.end()
+             && aIter != allRegionAlphaVolume.end();
+                ++vIter, ++aIter, ++numIter
+            )
             {
-                Info<< token::TAB << iter.key()
+                Info<< token::TAB << vIter.key()
-                    << token::TAB << iter().first()
+                    << token::TAB << vIter()
-                    << token::TAB << iter().second() << endl;
+                    << token::TAB << aIter()
+                    << token::TAB << numIter()
+                    << endl;
 
-                meshSumVol += iter().first();
+                meshSumVol += vIter();
-                alphaSumVol += iter().second();
+                alphaSumVol += aIter();
+                nCells += numIter();
             }
             Info<< token::TAB << "Total:"
                 << token::TAB << meshSumVol
-                << token::TAB << alphaSumVol << endl;
+                << token::TAB << alphaSumVol
+                << token::TAB << nCells
+                << endl;
             Info<< endl;
         }
 
 
 
-        // Mark all regions starting at patches
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-        // Count number of patch faces (just for initial sizing)
-        label nPatchFaces = 0;
-        forAll(patchNames_, i)
         {
-            const word& pName = patchNames_[i];
+            Info<< "Patch connected regions (liquid core):" << endl;
-            label patchI = mesh.boundaryMesh().findPatchID(pName);
+            Info<< token::TAB << "Region"
-            if (patchI == -1)
+                << token::TAB << "Volume(mesh)"
+                << token::TAB << "Volume(" << alpha.name() << "):"
+                << endl;
+            forAllConstIter(Map<label>, patchRegions, iter)
             {
-                WarningIn("regionSizeDistribution::write()")
+                label regionI = iter.key();
-                    << "Cannot find patch " << pName << ". Valid patches are "
+                Info<< token::TAB << iter.key()
-                    << mesh.boundaryMesh().names()
+                    << token::TAB << allRegionVolume[regionI]
-                    << endl;
+                    << token::TAB << allRegionAlphaVolume[regionI] << endl;
-            }
+
-            else
-            {
-                nPatchFaces += mesh.boundaryMesh()[patchI].size();
             }
+            Info<< endl;
         }
 
-        Map<label> keepRegions(nPatchFaces);
-        forAll(patchNames_, i)
         {
-            const word& pName = patchNames_[i];
+            Info<< "Background regions:" << endl;
+            Info<< token::TAB << "Region"
+                << token::TAB << "Volume(mesh)"
+                << token::TAB << "Volume(" << alpha.name() << "):"
+                << endl;
+            Map<scalar>::const_iterator vIter = allRegionVolume.begin();
+            Map<scalar>::const_iterator aIter = allRegionAlphaVolume.begin();
 
-            label patchI = mesh.boundaryMesh().findPatchID(pName);
+            for
-            if (patchI != -1)
-            {
-                const polyPatch& pp = mesh.boundaryMesh()[patchI];
-
-                // Collect all regions on the patch
-                const labelList& faceCells = pp.faceCells();
-
-                forAll(faceCells, i)
-                {
-                    keepRegions.insert
-                    (
-                        regions[faceCells[i]],
-                        Pstream::myProcNo()
-                    );
-                }
-            }
-        }
-
-
-        // Make sure all the processors have the same set of regions
-        Pstream::mapCombineGather(keepRegions, minEqOp<label>());
-        Pstream::mapCombineScatter(keepRegions);
-
-        Info<< "Patch connected regions (liquid core):" << endl;
-        forAllConstIter(Map<label>, keepRegions, iter)
-        {
-            label regionI = iter.key();
-            Pair<scalar>& vols = regionVolume[regionI];
-            Info<< token::TAB << iter.key()
-                << token::TAB << vols.first()
-                << token::TAB << vols.second() << endl;
-
-        }
-        Info<< endl;
-
-        Info<< "Background regions:" << endl;
-        forAllConstIter(Map<Pair<scalar> >, regionVolume, iter)
-        {
-            if
             (
-               !keepRegions.found(iter.key())
+                ;
-             && iter().first() >= volFraction_*meshVol
+                vIter != allRegionVolume.end()
+             && aIter != allRegionAlphaVolume.end();
+                ++vIter, ++aIter
             )
             {
-                Info<< token::TAB << iter.key()
+                if
-                    << token::TAB << iter().first()
+                (
-                    << token::TAB << iter().second() << endl;
+                   !patchRegions.found(vIter.key())
+                 && vIter() >= maxDropletVol
+                )
+                {
+                    Info<< token::TAB << vIter.key()
+                        << token::TAB << vIter()
+                        << token::TAB << aIter() << endl;
+                }
             }
+            Info<< endl;
         }
-        Info<< endl;
+
 
 
         // Split alpha field
@@ -400,185 +658,197 @@ void Foam::regionSizeDistribution::write()
         //  - liquidCore            : region connected to inlet patches
         //  - per region a volume   : for all other regions
         //  - backgroundAlpha       : remaining alpha
+        writeAlphaFields(regions, patchRegions, allRegionVolume, alpha);
 
 
-        // Construct field
+        // Extract droplet-only allRegionVolume, i.e. delete liquid core
-        volScalarField liquidCore
+        // (patchRegions) and background regions from maps.
-        (
+        // Note that we have to use mesh volume (allRegionVolume) and not
-            IOobject
+        // allRegionAlphaVolume since background might not have alpha in it.
-            (
+        forAllIter(Map<scalar>, allRegionVolume, vIter)
-                alphaName_ + "_liquidCore",
-                obr_.time().timeName(),
-                obr_,
-                IOobject::NO_READ
-            ),
-            alpha,
-            fvPatchField<scalar>::calculatedType()
-        );
-
-        volScalarField backgroundAlpha
-        (
-            IOobject
-            (
-                alphaName_ + "_background",
-                obr_.time().timeName(),
-                obr_,
-                IOobject::NO_READ
-            ),
-            alpha,
-            fvPatchField<scalar>::calculatedType()
-        );
-
-
-        // Knock out any cell not in keepRegions
-        forAll(liquidCore, cellI)
         {
-            label regionI = regions[cellI];
+            label regionI = vIter.key();
-            if (keepRegions.found(regionI))
+            if
+            (
+                patchRegions.found(regionI)
+             || vIter() >= maxDropletVol
+            )
             {
-                backgroundAlpha[cellI] = 0;
+                allRegionVolume.erase(vIter);
-            }
+                allRegionAlphaVolume.erase(regionI);
-            else
+                allRegionNumCells.erase(regionI);
-            {
-                liquidCore[cellI] = 0;
-
-                scalar regionVol = regionVolume[regionI].first();
-                if (regionVol < volFraction_*meshVol)
-                {
-                    backgroundAlpha[cellI] = 0;
-                }
             }
         }
-        liquidCore.correctBoundaryConditions();
-        backgroundAlpha.correctBoundaryConditions();
 
-        Info<< "Volume of liquid-core = "
+        if (allRegionVolume.size())
-            << fvc::domainIntegrate(liquidCore).value()
-            << endl;
-
-        Info<< "Writing liquid-core field to " << liquidCore.name() << endl;
-        liquidCore.write();
-
-        Info<< "Volume of background = "
-            << fvc::domainIntegrate(backgroundAlpha).value()
-            << endl;
-
-        Info<< "Writing background field to " << backgroundAlpha.name() << endl;
-        backgroundAlpha.write();
-
-
-
-        // Collect histogram
-        if (Pstream::master())
         {
-            DynamicList<scalar> diameters(regionVolume.size());
+            // Construct mids of bins for plotting
-            forAllConstIter(Map<Pair<scalar> >, regionVolume, iter)
+            pointField xBin(nBins_);
+
+            scalar x = 0.5*delta;
+            forAll(xBin, i)
             {
-                if (!keepRegions.found(iter.key()))
+                xBin[i] = point(x, 0, 0);
-                {
+                x += delta;
-                    if (iter().first() < volFraction_*meshVol)
-                    {
-                        scalar v = iter().second();
-                      //scalar diam = Foam::cbrt(v*6/mathematicalConstant::pi);
-                        scalar diam =
-                            Foam::cbrt(v*6/constant::mathematical::pi);
-                        diameters.append(diam);
-                    }
-                }
             }
 
-            if (diameters.size())
+            const coordSet coords("diameter", "x", xBin, mag(xBin));
+
+
+            // Get in region order the alpha*volume and diameter
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+            const labelList sortedRegions = allRegionAlphaVolume.sortedToc();
+
+            scalarField sortedVols
+            (
+                extractData
+                (
+                    sortedRegions,
+                    allRegionAlphaVolume
+                )
+            );
+
+            // Calculate the diameters
+            scalarField sortedDiameters(sortedVols.size());
+            forAll(sortedDiameters, i)
             {
-                scalar maxDiam = max(diameters);
+                sortedDiameters[i] = Foam::cbrt
-                scalar minDiam = 0.0;
-
-                Info<< "Maximum diameter:" << maxDiam << endl;
-
-                Histogram<List<scalar> > bins
                 (
-                    minDiam,
+                    sortedVols[i]
-                    maxDiam,
+                   *6/constant::mathematical::pi
-                    nBins_,
-                    diameters
                 );
+            }
 
-                /* 1.7.x
+            // Determine the bin index for all the diameters
-                scalarField xBin(nBins_);
+            labelList indices(sortedDiameters.size());
+            forAll(sortedDiameters, i)
+            {
+                indices[i] = (sortedDiameters[i]-minDiam_)/delta;
+            }
 
-                scalar dx = (maxDiam-minDiam)/nBins_;
+            // Calculate the counts per diameter bin
-                scalar x = 0.5*dx;
+            scalarField binCount(nBins_, 0.0);
-                forAll(bins.counts(), i)
+            forAll(sortedDiameters, i)
+            {
+                binCount[indices[i]] += 1.0;
+            }
+
+            // Write counts
+            if (Pstream::master())
+            {
+                writeGraph(coords, "count", binCount);
+            }
+
+            // Write to screen
+            {
+                Info<< "Bins:" << endl;
+                Info<< token::TAB << "Bin"
+                    << token::TAB << "Min diameter"
+                    << token::TAB << "Count:"
+                    << endl;
+
+                scalar diam = 0.0;
+                forAll(binCount, binI)
                 {
-                    xBin[i] = x;
+                    Info<< token::TAB << binI
-                    x += dx;
+                        << token::TAB << diam
+                        << token::TAB << binCount[binI] << endl;
+                    diam += delta;
                 }
+                Info<< endl;
+            }
 
-                scalarField normalisedCount(bins.counts().size());
+
-                forAll(bins.counts(), i)
+            // Write average and deviation of droplet volume.
+            writeGraphs
+            (
+                "volume",           // name of field
+                indices,            // per region the bin index
+                sortedVols,         // per region field data
+                binCount,           // per bin number of regions
+                coords              // graph data for bins
+            );
+
+            // Collect some more field
+            {
+                wordList scalarNames(obr_.names(volScalarField::typeName));
+                labelList selected = findStrings(fields_, scalarNames);
+
+                forAll(selected, i)
                 {
-                    normalisedCount[i] = 1.0*bins.counts()[i];
+                    const word& fldName = scalarNames[selected[i]];
+                    Info<< "Scalar field " << fldName << endl;
+
+                    const scalarField& fld = obr_.lookupObject
+                    <
+                        volScalarField
+                    >(fldName).internalField();
+
+                    writeGraphs
+                    (
+                        fldName,            // name of field
+                        alphaVol*fld,       // per cell field data
+
+                        regions,            // per cell the region(=droplet)
+                        sortedRegions,      // valid regions in sorted order
+                        1.0/sortedVols,     // per region normalisation
+
+                        indices,            // index of bin for each region
+                        binCount,           // per bin number of regions
+                        coords              // graph data for bins
+                    );
                 }
+            }
+            {
+                wordList vectorNames(obr_.names(volVectorField::typeName));
+                labelList selected = findStrings(fields_, vectorNames);
 
-                const coordSet coords
+                forAll(selected, i)
-                (
-                    "diameter",
-                    "x",
-                    xBin
-                );
-                */
-
-                pointField xBin(nBins_);
-                scalar dx = (maxDiam - minDiam)/nBins_;
-                scalar x = 0.5*dx;
-                forAll(bins.counts(), i)
                 {
-                    xBin[i] = point(x, 0, 0);
+                    const word& fldName = vectorNames[selected[i]];
-                    x += dx;
+                    Info<< "Vector field " << fldName << endl;
+
+                    const vectorField& fld = obr_.lookupObject
+                    <
+                        volVectorField
+                    >(fldName).internalField();
+
+
+                    // Components
+
+                    for (direction cmp = 0; cmp < vector::nComponents; cmp++)
+                    {
+                        writeGraphs
+                        (
+                            fldName + vector::componentNames[cmp],
+                            alphaVol*fld.component(cmp),// per cell field data
+
+                            regions,        // per cell the region(=droplet)
+                            sortedRegions,  // valid regions in sorted order
+                            1.0/sortedVols, // per region normalisation
+
+                            indices,        // index of bin for each region
+                            binCount,       // per bin number of regions
+                            coords          // graph data for bins
+                        );
+                    }
+
+                    // Magnitude
+                    writeGraphs
+                    (
+                        fldName + "mag",    // name of field
+                        alphaVol*mag(fld),  // per cell field data
+
+                        regions,            // per cell the region(=droplet)
+                        sortedRegions,      // valid regions in sorted order
+                        1.0/sortedVols,     // per region normalisation
+
+                        indices,            // index of bin for each region
+                        binCount,           // per bin number of regions
+                        coords              // graph data for bins
+                    );
                 }
-
-                scalarField normalisedCount(bins.counts().size());
-                forAll(bins.counts(), i)
-                {
-                    normalisedCount[i] = 1.0*bins.counts()[i];
-                }
-
-                const coordSet coords
-                (
-                    "diameter",
-                    "x",
-                    xBin,
-                    mag(xBin)
-                );
-                const wordList valNames(1, "count");
-
-
-                fileName outputPath;
-                if (Pstream::parRun())
-                {
-                    outputPath = mesh.time().path()/".."/name_;
-                }
-                else
-                {
-                    outputPath = mesh.time().path()/name_;
-                }
-
-                if (mesh.name() != fvMesh::defaultRegion)
-                {
-                    outputPath = outputPath/mesh.name();
-                }
-
-                mkDir(outputPath/mesh.time().timeName());
-                OFstream str
-                (
-                    outputPath
-                  / mesh.time().timeName()
-                  / formatterPtr_().getFileName(coords, valNames)
-                );
-                Info<< "Writing distribution to " << str.name() << endl;
-
-                List<const scalarField*> valPtrs(1);
-                valPtrs[0] = &normalisedCount;
-                formatterPtr_().write(coords, valNames, valPtrs, str);
             }
         }
     }

src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H

@@ -25,11 +25,67 @@ Class
     Foam::regionSizeDistribution
 
 Description
-    Looks up a field, interpolates it to the faces and determines a connected
+    Droplet size distribution calculation.
-    region from a patch where the field is above a certain value.
+
-    - Writes a field containing all regions starting at given patch
+    Looks up a void-fraction (alpha) field and splits the mesh into regions
-      ('liquid core')
+    based on where the field is below the threshold value. These
-    - All other regions are summed for volume and a histogram is calculated.
+    regions ("droplets") can now be analysed.
+
+    Regions:
+    - (debug) write regions as a volScalarField
+    - (debug) print for all regions the sum of volume and alpha*volume
+    - print the regions connected to a user-defined set of patches.
+      (in spray calculation these form the liquid core)
+    - print the regions with too large volume. These are the 'background'
+      regions.
+
+    Fields:
+    - write volScalarField alpha_liquidCore : alpha with outside liquid core
+                                              set to 0.
+                           alpha_background : alpha with outside background
+                                              set to 0.
+
+    Histogram:
+    - determine histogram of diameter (given minDiameter, maxDiameter, nBins)
+    - write graph of number of droplets per bin
+    - write graph of sum, average and deviation of droplet volume per bin
+    - write graph of sum, average and deviation of user-defined fields. For
+      volVectorFields these are those of the 3 components and the magnitude.
+
+    Sample input:
+
+    functions
+    {
+        regionSizeDistribution
+        {
+            type            regionSizeDistribution;
+
+            outputControl   timeStep;
+            outputInterval  1;
+
+            // Field to determine regions from
+            field           alpha;
+            // Patches that provide the liquid core
+            patches         (inlet);
+            // Delimit alpha regions
+            threshold       0.4;
+
+            // Fields to sample (no need to include alpha)
+            fields          (p U);
+
+            // Number of bins for histogram
+            nBins           100;
+            // Max droplet diameter
+            maxDiameter     0.5e-4;
+            //// Min droplet diameter (default is 0)
+            //minDiameter     0;
+
+            // Writing format
+            setFormat       gnuplot;
+        }
+    }
+
+
 
 SourceFiles
     regionSizeDistribution.C
@@ -41,6 +97,9 @@ SourceFiles
 
 #include "pointFieldFwd.H"
 #include "writer.H"
+#include "Map.H"
+#include "volFieldsFwd.H"
+#include "wordReList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,6 +110,8 @@ namespace Foam
 class objectRegistry;
 class dictionary;
 class mapPolyMesh;
+class regionSplit;
+class polyMesh;
 
 /*---------------------------------------------------------------------------*\
                            Class regionSizeDistribution Declaration
@@ -72,23 +133,85 @@ class regionSizeDistribution
         word alphaName_;
 
         //- Patches to walk from
-        wordList patchNames_;
+        wordReList patchNames_;
 
         //- Clip value
         scalar threshold_;
 
-        //- Background region volFraction
+        //- Maximum droplet diameter
-        scalar volFraction_;
+        scalar maxDiam_;
+
+        //- Minimum droplet diameter
+        scalar minDiam_;
 
         //- Mumber of bins
         label nBins_;
 
+        //- Names of fields to sample on regions
+        wordReList fields_;
+
         //- Output formatter to write
         autoPtr<writer<scalar> > formatterPtr_;
 
 
     // Private Member Functions
 
+        template<class Type>
+        Map<Type> regionSum(const regionSplit&, const Field<Type>&) const;
+
+        //- Get data in order
+        template<class Type>
+        List<Type> extractData(const UList<label>& keys, const Map<Type>&)
+        const;
+
+        void writeGraph
+        (
+            const coordSet& coords,
+            const word& valueName,
+            const scalarField& values
+        ) const;
+
+        //- Write volfields with the parts of alpha which are not
+        //  droplets (liquidCore, backGround)
+        void writeAlphaFields
+        (
+            const regionSplit& regions,
+            const Map<label>& keepRegions,
+            const Map<scalar>& regionVolume,
+            const volScalarField& alpha
+        ) const;
+
+        //- Mark all regions starting at patches
+        Map<label> findPatchRegions(const polyMesh&, const regionSplit&) const;
+
+        //- Helper: divide if denom != 0
+        static tmp<scalarField> divide(const scalarField&, const scalarField&);
+
+        //- Given per-region data calculate per-bin average/deviation and graph
+        void writeGraphs
+        (
+            const word& fieldName,              // name of field
+            const labelList& indices,           // index of bin for each region
+            const scalarField& sortedField,     // per region field data
+            const scalarField& binCount,        // per bin number of regions
+            const coordSet& coords              // graph data for bins
+        ) const;
+
+        //- Given per-cell data calculate per-bin average/deviation and graph
+        void writeGraphs
+        (
+            const word& fieldName,              // name of field
+            const scalarField& cellField,       // per cell field data
+
+            const regionSplit& regions,         // per cell the region(=droplet)
+            const labelList& sortedRegions,     // valid regions in sorted order
+            const scalarField& sortedNormalisation,
+
+            const labelList& indices,           // index of bin for each region
+            const scalarField& binCount,        // per bin number of regions
+            const coordSet& coords              // graph data for bins
+        ) const;
+
         //- Disallow default bitwise copy construct
         regionSizeDistribution(const regionSizeDistribution&);
 
@@ -156,6 +279,12 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef NoRepository
+#   include "regionSizeDistributionTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //

src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistributionTemplates.C

@@ -0,0 +1,81 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "regionSizeDistribution.H"
+#include "regionSplit.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::Map<Type> Foam::regionSizeDistribution::regionSum
+(
+    const regionSplit& regions,
+    const Field<Type>& fld
+) const
+{
+    // Per region the sum of fld
+    Map<Type> regionToSum(regions.nRegions()/Pstream::nProcs());
+
+    forAll(fld, cellI)
+    {
+        label regionI = regions[cellI];
+
+        typename Map<Type>::iterator fnd = regionToSum.find(regionI);
+        if (fnd == regionToSum.end())
+        {
+            regionToSum.insert(regionI, fld[cellI]);
+        }
+        else
+        {
+            fnd() += fld[cellI];
+        }
+    }
+    Pstream::mapCombineGather(regionToSum, plusEqOp<Type>());
+    Pstream::mapCombineScatter(regionToSum);
+
+    return regionToSum;
+}
+
+
+// Get data in sortedToc order
+template<class Type>
+Foam::List<Type> Foam::regionSizeDistribution::extractData
+(
+    const UList<label>& keys,
+    const Map<Type>& regionData
+) const
+{
+    List<Type> sortedData(keys.size());
+
+    forAll(keys, i)
+    {
+        sortedData[i] = regionData[keys[i]];
+    }
+    return sortedData;
+}
+
+
+// ************************************************************************* //

src/postProcessing/functionObjects/forces/forces/forces.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -66,14 +66,14 @@ Foam::tmp<Foam::volSymmTensorField> Foam::forces::devRhoReff() const
         const compressible::LESModel& les =
         obr_.lookupObject<compressible::LESModel>("LESProperties");
 
-        return les.devRhoBeff();
+        return les.devRhoReff();
     }
     else if (obr_.foundObject<incompressible::LESModel>("LESProperties"))
     {
         const incompressible::LESModel& les
             = obr_.lookupObject<incompressible::LESModel>("LESProperties");
 
-        return rho()*les.devBeff();
+        return rho()*les.devReff();
     }
     else if (obr_.foundObject<basicThermo>("thermophysicalProperties"))
     {

src/sampling/sampledSurface/writers/vtk/vtkSurfaceWriter.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -170,6 +170,7 @@ namespace Foam
         {
             const symmTensor& v = values[elemI];
             os  << float(v[0]) << ' ' << float(v[1]) << ' ' << float(v[2])
+                << ' '
                 << float(v[3]) << ' ' << float(v[4]) << ' ' << float(v[5])
                 << nl;
 
@@ -190,7 +191,9 @@ namespace Foam
         {
             const tensor& v = values[elemI];
             os  << float(v[0]) << ' ' << float(v[1]) << ' ' << float(v[2])
+                << ' '
                 << float(v[3]) << ' ' << float(v[4]) << ' ' << float(v[5])
+                << ' '
                 << float(v[6]) << ' ' << float(v[7]) << ' ' << float(v[8])
                 << nl;
         }

src/thermophysicalModels/radiationModels/derivedFvPatchFields/greyDiffusiveViewFactor/greyDiffusiveViewFactorFixedValueFvPatchScalarField.C

@@ -128,10 +128,14 @@ greyDiffusiveViewFactorFixedValueFvPatchScalarField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-
 void Foam::radiation::greyDiffusiveViewFactorFixedValueFvPatchScalarField::
 updateCoeffs()
 {
+    if (this->updated())
+    {
+        return;
+    }
+
     // Do nothing
 
     if (debug)
@@ -149,6 +153,7 @@ updateCoeffs()
             << endl;
     }
 
+    fixedValueFvPatchScalarField::updateCoeffs();
 }
 
 

src/turbulenceModels/compressible/RAS/derivedFvPatchFields/convectiveHeatTransfer/convectiveHeatTransferFvPatchScalarField.C

@@ -142,6 +142,8 @@ void convectiveHeatTransferFvPatchScalarField::updateCoeffs()
             htc[faceI] = 0.037*pow(Re, 0.8)*cbrt(Pr[faceI])*kappaw[faceI]/L_;
         }
     }
+
+    fixedValueFvPatchScalarField::updateCoeffs();
 }
 
 

src/turbulenceModels/incompressible/LES/GenEddyVisc/GenEddyVisc.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -87,17 +87,34 @@ tmp<volSymmTensorField> GenEddyVisc::B() const
 }
 
 
-tmp<volSymmTensorField> GenEddyVisc::devBeff() const
+tmp<volSymmTensorField> GenEddyVisc::devReff() const
 {
     return -nuEff()*dev(twoSymm(fvc::grad(U())));
 }
 
 
-tmp<fvVectorMatrix> GenEddyVisc::divDevBeff(volVectorField& U) const
+tmp<fvVectorMatrix> GenEddyVisc::divDevReff(volVectorField& U) const
 {
     return
     (
-      - fvm::laplacian(nuEff(), U) - fvc::div(nuEff()*dev(T(fvc::grad(U))))
+      - fvm::laplacian(nuEff(), U)
+      - fvc::div(nuEff()*dev(T(fvc::grad(U))))
+    );
+}
+
+
+tmp<fvVectorMatrix> GenEddyVisc::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
     );
 }
 

src/turbulenceModels/incompressible/LES/GenEddyVisc/GenEddyVisc.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -118,11 +118,19 @@ public:
 
         //- Return the effective sub-grid turbulence stress tensor
         //  including the laminar stress
-        virtual tmp<volSymmTensorField> devBeff() const;
+        virtual tmp<volSymmTensorField> devReff() const;
 
         //- Return the deviatoric part of the effective sub-grid
         //  turbulence stress tensor including the laminar stress
-        virtual tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const;
+        virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
+
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
 
         //- Correct Eddy-Viscosity and related properties
         virtual void correct(const tmp<volTensorField>& gradU);

src/turbulenceModels/incompressible/LES/GenSGSStress/GenSGSStress.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -115,7 +115,7 @@ GenSGSStress::GenSGSStress
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-tmp<volSymmTensorField> GenSGSStress::devBeff() const
+tmp<volSymmTensorField> GenSGSStress::devReff() const
 {
     return tmp<volSymmTensorField>
     (
@@ -135,7 +135,7 @@ tmp<volSymmTensorField> GenSGSStress::devBeff() const
 }
 
 
-tmp<fvVectorMatrix> GenSGSStress::divDevBeff
+tmp<fvVectorMatrix> GenSGSStress::divDevReff
 (
     volVectorField& U
 ) const
@@ -164,6 +164,38 @@ tmp<fvVectorMatrix> GenSGSStress::divDevBeff
 }
 
 
+tmp<fvVectorMatrix> GenSGSStress::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    if (couplingFactor_.value() > 0.0)
+    {
+        return
+        (
+            fvc::div(rho*B_ + couplingFactor_*rho*nuSgs_*fvc::grad(U))
+          + fvc::laplacian
+            (
+                (1.0 - couplingFactor_)*rho*nuSgs_, U, "laplacian(muEff,U)"
+            )
+          - fvm::laplacian(muEff, U)
+        );
+    }
+    else
+    {
+        return
+        (
+            fvc::div(rho*B_)
+          + fvc::laplacian(rho*nuSgs_, U, "laplacian(muEff,U)")
+          - fvm::laplacian(muEff, U)
+        );
+    }
+}
+
+
 bool GenSGSStress::read()
 {
     if (LESModel::read())

src/turbulenceModels/incompressible/LES/GenSGSStress/GenSGSStress.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -128,11 +128,19 @@ public:
 
         //- Return the effective sub-grid turbulence stress tensor
         //  including the laminar stress
-        virtual tmp<volSymmTensorField> devBeff() const;
+        virtual tmp<volSymmTensorField> devReff() const;
 
-        //- Returns div(B).
+        //- Return the deviatoric part of the effective sub-grid
-        // This is the additional term due to the filtering of the NSE.
+        //  turbulence stress tensor including the laminar stress
-        virtual tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const;
+        virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
+
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
 
         //- Read LESProperties dictionary
         virtual bool read();

src/turbulenceModels/incompressible/LES/LESModel/LESModel.H

@@ -202,14 +202,6 @@ public:
         //- Return the sub-grid stress tensor.
         virtual tmp<volSymmTensorField> B() const = 0;
 
-        //- Return the deviatoric part of the effective sub-grid
-        //  turbulence stress tensor including the laminar stress
-        virtual tmp<volSymmTensorField> devBeff() const = 0;
-
-        //- Returns div(dev(Beff)).
-        //  This is the additional term due to the filtering of the NSE.
-        virtual tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const = 0;
-
 
         // RAS compatibility functions for the turbulenceModel base class
 
@@ -225,18 +217,6 @@ public:
                 return B();
             }
 
-            //- Return the effective stress tensor including the laminar stress
-            virtual tmp<volSymmTensorField> devReff() const
-            {
-                return devBeff();
-            }
-
-            //- Return the source term for the momentum equation
-            virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const
-            {
-                return divDevBeff(U);
-            }
-
 
         //- Correct Eddy-Viscosity and related properties.
         //  This calls correct(const tmp<volTensorField>& gradU) by supplying

src/turbulenceModels/incompressible/LES/Smagorinsky2/Smagorinsky2.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -81,7 +81,7 @@ tmp<volSymmTensorField> Smagorinsky2::B() const
 }
 
 
-tmp<fvVectorMatrix> Smagorinsky2::divDevBeff
+tmp<fvVectorMatrix> Smagorinsky2::divDevReff
 (
     volVectorField& U
 ) const
@@ -101,6 +101,28 @@ tmp<fvVectorMatrix> Smagorinsky2::divDevBeff
 }
 
 
+tmp<fvVectorMatrix> Smagorinsky2::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volTensorField gradU(fvc::grad(U));
+
+    volSymmTensorField aniMuEff
+    (
+        "muEff",
+        I*(rho*nuEff()) + (cD2_*rho*delta())*symm(gradU)
+    );
+
+    return
+    (
+      - fvm::laplacian(aniMuEff, U)
+      - fvc::div(rho*nuEff()*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool Smagorinsky2::read()
 {
     if (Smagorinsky::read())

src/turbulenceModels/incompressible/LES/Smagorinsky2/Smagorinsky2.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -109,9 +109,17 @@ public:
         //- Return B.
         virtual tmp<volSymmTensorField> B() const;
 
-        //- Returns div(B).
+        //- Return the deviatoric part of the effective sub-grid
-        // This is the additional term due to the filtering of the NSE.
+        //  turbulence stress tensor including the laminar stress
-        virtual tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const;
+        virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
+
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
 
         //- Read LESProperties dictionary
         virtual bool read();

src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -338,17 +338,34 @@ tmp<volSymmTensorField> SpalartAllmaras::B() const
 }
 
 
-tmp<volSymmTensorField> SpalartAllmaras::devBeff() const
+tmp<volSymmTensorField> SpalartAllmaras::devReff() const
 {
     return -nuEff()*dev(twoSymm(fvc::grad(U())));
 }
 
 
-tmp<fvVectorMatrix> SpalartAllmaras::divDevBeff(volVectorField& U) const
+tmp<fvVectorMatrix> SpalartAllmaras::divDevReff(volVectorField& U) const
 {
     return
     (
-      - fvm::laplacian(nuEff(), U) - fvc::div(nuEff()*dev(T(fvc::grad(U))))
+      - fvm::laplacian(nuEff(), U)
+      - fvc::div(nuEff()*dev(T(fvc::grad(U))))
+    );
+}
+
+
+tmp<fvVectorMatrix> SpalartAllmaras::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
     );
 }
 

src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -172,11 +172,19 @@ public:
 
         //- Return the effective sub-grid turbulence stress tensor
         //  including the laminar stress
-        virtual tmp<volSymmTensorField> devBeff() const;
+        virtual tmp<volSymmTensorField> devReff() const;
 
-        //- Return the deviatoric part of the divergence of Beff
+        //- Return the deviatoric part of the effective sub-grid
-        //  i.e. the additional term in the filtered NSE.
+        //  turbulence stress tensor including the laminar stress
-        virtual tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const;
+        virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
+
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
 
         //- Correct nuTilda and related properties
         virtual void correct(const tmp<volTensorField>& gradU);

src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -388,7 +388,7 @@ void kOmegaSSTSAS::correct(const tmp<volTensorField>& gradU)
           - fvm::Sp(fvc::div(phi()), omega_)
           - fvm::laplacian(DomegaEff(F1), omega_)
         ==
-            gamma(F1)*0.5*S2
+            gamma(F1)*S2
           - fvm::Sp(beta(F1)*omega_, omega_)
           - fvm::SuSp       // cross diffusion term
             (
@@ -425,17 +425,34 @@ tmp<volSymmTensorField> kOmegaSSTSAS::B() const
 }
 
 
-tmp<volSymmTensorField> kOmegaSSTSAS::devBeff() const
+tmp<volSymmTensorField> kOmegaSSTSAS::devReff() const
 {
     return -nuEff()*dev(twoSymm(fvc::grad(U())));
 }
 
 
-tmp<fvVectorMatrix> kOmegaSSTSAS::divDevBeff(volVectorField& U) const
+tmp<fvVectorMatrix> kOmegaSSTSAS::divDevReff(volVectorField& U) const
 {
     return
     (
-      - fvm::laplacian(nuEff(), U) - fvc::div(nuEff()*dev(T(fvc::grad(U))))
+      - fvm::laplacian(nuEff(), U)
+      - fvc::div(nuEff()*dev(T(fvc::grad(U))))
+    );
+}
+
+
+tmp<fvVectorMatrix> kOmegaSSTSAS::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
     );
 }
 

src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -249,11 +249,19 @@ public:
 
         //- Return the effective sub-grid turbulence stress tensor
         //  including the laminar stress
-        virtual tmp<volSymmTensorField> devBeff() const;
+        virtual tmp<volSymmTensorField> devReff() const;
 
-        //- Return the deviatoric part of the divergence of Beff
+        //- Return the deviatoric part of the effective sub-grid
-        //  i.e. the additional term in the filtered NSE.
+        //  turbulence stress tensor including the laminar stress
-        virtual tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const;
+        virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
+
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
 
         //- Solve the turbulence equations (k-w) and correct the turbulence
         //  viscosity

src/turbulenceModels/incompressible/LES/laminar/laminar.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -136,17 +136,34 @@ tmp<volSymmTensorField> laminar::B() const
 }
 
 
-tmp<volSymmTensorField> laminar::devBeff() const
+tmp<volSymmTensorField> laminar::devReff() const
 {
     return -nu()*dev(twoSymm(fvc::grad(U())));
 }
 
 
-tmp<fvVectorMatrix> laminar::divDevBeff(volVectorField& U) const
+tmp<fvVectorMatrix> laminar::divDevReff(volVectorField& U) const
 {
     return
     (
-      - fvm::laplacian(nu(), U) - fvc::div(nu()*dev(T(fvc::grad(U))))
+      - fvm::laplacian(nu(), U)
+      - fvc::div(nu()*dev(T(fvc::grad(U))))
+    );
+}
+
+
+tmp<fvVectorMatrix> laminar::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
     );
 }
 

src/turbulenceModels/incompressible/LES/laminar/laminar.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -104,13 +104,21 @@ public:
         //- Return the sub-grid stress tensor B.
         virtual tmp<volSymmTensorField> B() const;
 
+        //- Return the effective sub-grid turbulence stress tensor
+        //  including the laminar stress
+        virtual tmp<volSymmTensorField> devReff() const;
+
         //- Return the deviatoric part of the effective sub-grid
         //  turbulence stress tensor including the laminar stress
-        virtual tmp<volSymmTensorField> devBeff() const;
+        virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
-        //- Return the deviatoric part of the divergence of Beff
+        //- Return the deviatoric part of the effective sub-grid
-        //  i.e. the additional term in the filtered NSE.
+        //  turbulence stress tensor including the laminar stress
-        virtual tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const;
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
 
         //- Read LESProperties dictionary
         bool read();

src/turbulenceModels/incompressible/LES/mixedSmagorinsky/mixedSmagorinsky.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -91,25 +91,39 @@ tmp<volSymmTensorField> mixedSmagorinsky::B() const
 }
 
 
-tmp<volSymmTensorField> mixedSmagorinsky::devBeff() const
+tmp<volSymmTensorField> mixedSmagorinsky::devReff() const
 {
     return
     (
-        scaleSimilarity::devBeff()
+        scaleSimilarity::devReff()
-      + Smagorinsky::devBeff()
+      + Smagorinsky::devReff()
     );
 }
 
 
-tmp<fvVectorMatrix> mixedSmagorinsky::divDevBeff
+tmp<fvVectorMatrix> mixedSmagorinsky::divDevReff
 (
     volVectorField& U
 ) const
 {
     return
     (
-        scaleSimilarity::divDevBeff(U)
+        scaleSimilarity::divDevReff(U)
-      + Smagorinsky::divDevBeff(U)
+      + Smagorinsky::divDevReff(U)
+    );
+}
+
+
+tmp<fvVectorMatrix> mixedSmagorinsky::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    return
+    (
+        scaleSimilarity::divDevRhoReff(rho, U)
+      + Smagorinsky::divDevRhoReff(rho, U)
     );
 }
 

src/turbulenceModels/incompressible/LES/mixedSmagorinsky/mixedSmagorinsky.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -129,11 +129,19 @@ public:
 
         //- Return the effective sub-grid turbulence stress tensor
         //  including the laminar stress
-        virtual tmp<volSymmTensorField> devBeff() const;
+        virtual tmp<volSymmTensorField> devReff() const;
 
-        //- Implementation of div(B). This is necessary to override
+        //- Return the deviatoric part of the effective sub-grid
-        // (and include) the div(B) terms from both the parent classes.
+        //  turbulence stress tensor including the laminar stress
-        virtual tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const;
+        virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
+
+        //- Return the deviatoric part of the effective sub-grid
+        //  turbulence stress tensor including the laminar stress
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
 
         //- Correct Eddy-Viscosity and related properties
         virtual void correct(const tmp<volTensorField>& gradU);

src/turbulenceModels/incompressible/LES/scaleSimilarity/scaleSimilarity.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -79,15 +79,25 @@ tmp<volSymmTensorField> scaleSimilarity::B() const
 }
 
 
-tmp<volSymmTensorField> scaleSimilarity::devBeff() const
+tmp<volSymmTensorField> scaleSimilarity::devReff() const
 {
     return dev(B());
 }
 
 
-tmp<fvVectorMatrix> scaleSimilarity::divDevBeff(volVectorField& U) const
+tmp<fvVectorMatrix> scaleSimilarity::divDevReff(volVectorField& U) const
 {
-    return fvm::Su(fvc::div(devBeff()), U);
+    return fvm::Su(fvc::div(devReff()), U);
+}
+
+
+tmp<fvVectorMatrix> scaleSimilarity::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    return fvm::Su(fvc::div(rho*devReff()), U);
 }
 
 

src/turbulenceModels/incompressible/LES/scaleSimilarity/scaleSimilarity.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -106,13 +106,21 @@ public:
         //- Return the sub-grid stress tensor.
         virtual tmp<volSymmTensorField> B() const;
 
+        //- Return the effective sub-grid turbulence stress tensor
+        //  including the laminar stress
+        virtual tmp<volSymmTensorField> devReff() const;
+
         //- Return the deviatoric part of the effective sub-grid
         //  turbulence stress tensor including the laminar stress
-        virtual tmp<volSymmTensorField> devBeff() const;
+        virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
-        //- Return the deviatoric part of the divergence of Beff
+        //- Return the deviatoric part of the effective sub-grid
-        //  i.e. the additional term in the filtered NSE.
+        //  turbulence stress tensor including the laminar stress
-        virtual tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const;
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
 
         //- Correct Eddy-Viscosity and related properties
         virtual void correct(const tmp<volTensorField>&);

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -259,6 +259,44 @@ tmp<fvVectorMatrix> LRR::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> LRR::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    if (couplingFactor_.value() > 0.0)
+    {
+        return
+        (
+            fvc::div
+            (
+                rho*R_ + couplingFactor_*(rho*nut_)*fvc::grad(U),
+                "div((rho*R))"
+            )
+          + fvc::laplacian
+            (
+                (1.0 - couplingFactor_)*rho*nut_,
+                U,
+                "laplacian(muEff,U)"
+            )
+          - fvm::laplacian(muEff, U)
+        );
+    }
+    else
+    {
+        return
+        (
+            fvc::div(rho*R_)
+          + fvc::laplacian(rho*nut_, U, "laplacian(muEff,U)")
+          - fvm::laplacian(muEff, U)
+        );
+    }
+}
+
+
 bool LRR::read()
 {
     if (RASModel::read())

src/turbulenceModels/incompressible/RAS/LRR/LRR.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -173,6 +173,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 

src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -204,6 +204,22 @@ tmp<fvVectorMatrix> LamBremhorstKE::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> LamBremhorstKE::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool LamBremhorstKE::read()
 {
     if (RASModel::read())

src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -147,6 +147,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -276,7 +276,12 @@ tmp<fvVectorMatrix> LaunderGibsonRSTM::divDevReff(volVectorField& U) const
         return
         (
             fvc::div(R_ + couplingFactor_*nut_*fvc::grad(U), "div(R)")
-          + fvc::laplacian((1.0-couplingFactor_)*nut_, U, "laplacian(nuEff,U)")
+          + fvc::laplacian
+            (
+                (1.0 - couplingFactor_)*nut_,
+                U,
+                "laplacian(nuEff,U)"
+            )
           - fvm::laplacian(nuEff(), U)
         );
     }
@@ -292,6 +297,44 @@ tmp<fvVectorMatrix> LaunderGibsonRSTM::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> LaunderGibsonRSTM::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    if (couplingFactor_.value() > 0.0)
+    {
+        return
+        (
+            fvc::div
+            (
+                rho*R_ + couplingFactor_*(rho*nut_)*fvc::grad(U),
+                "div((rho*R))"
+            )
+          + fvc::laplacian
+            (
+                (1.0 - couplingFactor_)*rho*nut_,
+                U,
+                "laplacian(muEff,U)"
+            )
+          - fvm::laplacian(muEff, U)
+        );
+    }
+    else
+    {
+        return
+        (
+            fvc::div(rho*R_)
+          + fvc::laplacian(rho*nut_, U, "laplacian(muEff,U)")
+          - fvm::laplacian(muEff, U)
+        );
+    }
+}
+
+
 bool LaunderGibsonRSTM::read()
 {
     if (RASModel::read())

src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -184,6 +184,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 

src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -210,6 +210,22 @@ tmp<fvVectorMatrix> LaunderSharmaKE::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> LaunderSharmaKE::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool LaunderSharmaKE::read()
 {
     if (RASModel::read())

src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -164,6 +164,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 

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

@@ -298,6 +298,23 @@ tmp<fvVectorMatrix> LienCubicKE::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> LienCubicKE::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+        fvc::div(rho*nonlinearStress_)
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool LienCubicKE::read()
 {
     if (RASModel::read())

src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.H

@@ -159,6 +159,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 

src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C

@@ -365,6 +365,23 @@ tmp<fvVectorMatrix> LienCubicKELowRe::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> LienCubicKELowRe::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+        fvc::div(rho*nonlinearStress_)
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool LienCubicKELowRe::read()
 {
     if (RASModel::read())

src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.H

@@ -186,6 +186,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 

src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -243,12 +243,27 @@ tmp<fvVectorMatrix> LienLeschzinerLowRe::divDevReff(volVectorField& U) const
     return
     (
       - fvm::laplacian(nuEff(), U)
-    //- (fvc::grad(U) & fvc::grad(nuEff()))
       - fvc::div(nuEff()*T(fvc::grad(U)))
     );
 }
 
 
+tmp<fvVectorMatrix> LienLeschzinerLowRe::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool LienLeschzinerLowRe::read()
 {
     if (RASModel::read())

src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -156,6 +156,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 

src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.C

@@ -287,6 +287,23 @@ tmp<fvVectorMatrix> NonlinearKEShih::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> NonlinearKEShih::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+        fvc::div(rho*nonlinearStress_)
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool NonlinearKEShih::read()
 {
     if (RASModel::read())

src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.H

@@ -162,6 +162,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 

src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -220,6 +220,22 @@ tmp<fvVectorMatrix> RNGkEpsilon::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> RNGkEpsilon::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool RNGkEpsilon::read()
 {
     if (RASModel::read())

src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -161,6 +161,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 

src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.C

@@ -337,6 +337,22 @@ tmp<fvVectorMatrix> SpalartAllmaras::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> SpalartAllmaras::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool SpalartAllmaras::read()
 {
     if (RASModel::read())

src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -184,6 +184,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.C

@@ -116,7 +116,8 @@ atmBoundaryLayerInletEpsilonFvPatchScalarField
 
     z_ /= mag(z_);
 
-    evaluate();
+    const vectorField& c = patch().Cf();
+    scalarField::operator=(pow3(Ustar_)/(kappa_*((c & z_) - zGround_ + z0_)));
 }
 
 
@@ -169,14 +170,6 @@ void atmBoundaryLayerInletEpsilonFvPatchScalarField::rmap
 }
 
 
-void atmBoundaryLayerInletEpsilonFvPatchScalarField::updateCoeffs()
-{
-    const vectorField& c = patch().Cf();
-    tmp<scalarField> coord = (c & z_);
-    scalarField::operator=(pow3(Ustar_)/(kappa_*(coord - zGround_ + z0_)));
-}
-
-
 void atmBoundaryLayerInletEpsilonFvPatchScalarField::write(Ostream& os) const
 {
     fvPatchScalarField::write(os);

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.H

@@ -176,7 +176,7 @@ public:
         // Access
 
             //- Return max value
-            scalarField Ustar() const
+            const scalarField& Ustar() const
             {
                 return Ustar_;
             }
@@ -204,12 +204,6 @@ public:
             );
 
 
-        // Evaluation functions
-
-            //- Update coefficients
-            virtual void updateCoeffs();
-
-
         //- Write
         virtual void write(Ostream&) const;
 };

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.C

@@ -120,7 +120,25 @@ atmBoundaryLayerInletVelocityFvPatchVectorField
         Ustar_[i] = kappa_*Uref_/(log((Href_  + z0_[i])/max(z0_[i] , 0.001)));
     }
 
-    evaluate();
+    const vectorField& c = patch().Cf();
+    const scalarField coord(c & z_);
+    scalarField Un(coord.size());
+
+    forAll(coord, i)
+    {
+        if ((coord[i] - zGround_[i]) < Href_)
+        {
+            Un[i] =
+                (Ustar_[i]/kappa_)
+              * log((coord[i] - zGround_[i] + z0_[i])/max(z0_[i], 0.001));
+        }
+        else
+        {
+            Un[i] = Uref_;
+        }
+    }
+
+    vectorField::operator=(n_*Un);
 }
 
 
@@ -174,32 +192,6 @@ void atmBoundaryLayerInletVelocityFvPatchVectorField::rmap
 }
 
 
-void atmBoundaryLayerInletVelocityFvPatchVectorField::updateCoeffs()
-{
-    const vectorField& c = patch().Cf();
-    const scalarField coord(c & z_);
-    scalarField Un(coord.size());
-
-    forAll(coord, i)
-    {
-        if ((coord[i] - zGround_[i]) < Href_)
-        {
-            Un[i] =
-                (Ustar_[i]/kappa_)
-              * log((coord[i] - zGround_[i] + z0_[i])/max(z0_[i], 0.001));
-        }
-        else
-        {
-            Un[i] = Uref_;
-        }
-    }
-
-    vectorField::operator=(n_*Un);
-
-    fixedValueFvPatchVectorField::updateCoeffs();
-}
-
-
 void atmBoundaryLayerInletVelocityFvPatchVectorField::write(Ostream& os) const
 {
     fvPatchVectorField::write(os);

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.H

@@ -180,23 +180,26 @@ public:
 
     // Member functions
 
-        //- Return Ustar
+        // Access
-        scalarField& Ustar()
-        {
-            return Ustar_;
-        }
 
-        //- Return flow direction
+            //- Return Ustar
-        vector& n()
+            const scalarField& Ustar() const
-        {
+            {
-            return n_;
+                return Ustar_;
-        }
+            }
+
+            //- Return flow direction
+            const vector& n() const
+            {
+                return n_;
+            }
+
+            //- Return z direction
+            const vector& z() const
+            {
+                return z_;
+            }
 
-        //- Return z direction
-        vector& z()
-        {
-            return z_;
-        }
 
         // Mapping functions
 
@@ -214,12 +217,6 @@ public:
             );
 
 
-        // Evaluation functions
-
-            //- Update coefficients
-            virtual void updateCoeffs();
-
-
         //- Write
         virtual void write(Ostream&) const;
 };

src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -192,6 +192,22 @@ tmp<fvVectorMatrix> kEpsilon::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> kEpsilon::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool kEpsilon::read()
 {
     if (RASModel::read())

src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -155,6 +155,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 

src/turbulenceModels/incompressible/RAS/kOmega/kOmega.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -201,6 +201,22 @@ tmp<fvVectorMatrix> kOmega::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> kOmega::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool kOmega::read()
 {
     if (RASModel::read())

src/turbulenceModels/incompressible/RAS/kOmega/kOmega.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -188,6 +188,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 

src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -308,6 +308,22 @@ tmp<fvVectorMatrix> kOmegaSST::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> kOmegaSST::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool kOmegaSST::read()
 {
     if (RASModel::read())

src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -258,6 +258,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 

src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.C

@@ -570,6 +570,22 @@ tmp<fvVectorMatrix> kkLOmega::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> kkLOmega::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool kkLOmega::read()
 {
     if (RASModel::read())

src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -276,6 +276,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 

src/turbulenceModels/incompressible/RAS/laminar/laminar.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -177,6 +177,22 @@ tmp<fvVectorMatrix> laminar::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> laminar::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool laminar::read()
 {
     return RASModel::read();

src/turbulenceModels/incompressible/RAS/laminar/laminar.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -105,6 +105,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Correct the laminar viscosity
         virtual void correct();
 

src/turbulenceModels/incompressible/RAS/qZeta/qZeta.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -262,6 +262,22 @@ tmp<fvVectorMatrix> qZeta::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> qZeta::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool qZeta::read()
 {
     if (RASModel::read())

src/turbulenceModels/incompressible/RAS/qZeta/qZeta.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -211,6 +211,13 @@ public:
         //- Return the source term for the momentum equation
         virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
 
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
 

src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -246,6 +246,22 @@ tmp<fvVectorMatrix> realizableKE::divDevReff(volVectorField& U) const
 }
 
 
+tmp<fvVectorMatrix> realizableKE::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    volScalarField muEff("muEff", rho*nuEff());
+
+    return
+    (
+      - fvm::laplacian(muEff, U)
+      - fvc::div(muEff*dev(T(fvc::grad(U))))
+    );
+}
+
+
 bool realizableKE::read()
 {
     if (RASModel::read())