From f7040f7cc1a0a0018b27e3080338c3777199de00 Mon Sep 17 00:00:00 2001
From: andy <a.heather@opencfd.co.uk>
Date: Thu, 29 Jan 2009 13:41:08 +0000
Subject: [PATCH 1/5] using DpDt and adding run-time modifiable timestep
 support

---
 .../solvers/heatTransfer/buoyantFoam/buoyantFoam.C         | 7 ++++++-
 .../solvers/heatTransfer/buoyantFoam/createFields.H        | 5 +++--
 applications/solvers/heatTransfer/buoyantFoam/hEqn.H       | 4 +---
 applications/solvers/heatTransfer/buoyantFoam/pEqn.H       | 2 +-
 4 files changed, 11 insertions(+), 7 deletions(-)

diff --git a/applications/solvers/heatTransfer/buoyantFoam/buoyantFoam.C b/applications/solvers/heatTransfer/buoyantFoam/buoyantFoam.C
index cafb677900..24ad89e4e3 100644
--- a/applications/solvers/heatTransfer/buoyantFoam/buoyantFoam.C
+++ b/applications/solvers/heatTransfer/buoyantFoam/buoyantFoam.C
@@ -48,6 +48,9 @@ int main(int argc, char *argv[])
 #   include "readEnvironmentalProperties.H"
 #   include "createFields.H"
 #   include "initContinuityErrs.H"
+#   include "readTimeControls.H"
+#   include "compressibleCourantNo.H"
+#   include "setInitialDeltaT.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -55,11 +58,13 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
+#       include "readTimeControls.H"
 #       include "readPISOControls.H"
 #       include "compressibleCourantNo.H"
-//#       include "setDeltaT.H"
+#       include "setDeltaT.H"
 
         runTime++;
+
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
 #       include "rhoEqn.H"
diff --git a/applications/solvers/heatTransfer/buoyantFoam/createFields.H b/applications/solvers/heatTransfer/buoyantFoam/createFields.H
index 3d604d4daf..9535718fbb 100644
--- a/applications/solvers/heatTransfer/buoyantFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantFoam/createFields.H
@@ -52,8 +52,9 @@
         )
     );
 
-    Info<< "Creating field dpdt\n" << endl;
-    volScalarField dpdt = fvc::ddt(p);
+    Info<< "Creating field DpDt\n" << endl;
+    volScalarField DpDt =
+        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
 
     Info<< "Calculating field g.h\n" << endl;
     volScalarField gh("gh", g & mesh.C());
diff --git a/applications/solvers/heatTransfer/buoyantFoam/hEqn.H b/applications/solvers/heatTransfer/buoyantFoam/hEqn.H
index 008e3b0f56..f1a87c6a79 100644
--- a/applications/solvers/heatTransfer/buoyantFoam/hEqn.H
+++ b/applications/solvers/heatTransfer/buoyantFoam/hEqn.H
@@ -5,9 +5,7 @@
       + fvm::div(phi, h)
       - fvm::laplacian(turbulence->alphaEff(), h)
      ==
-        dpdt
-      + fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p))
-      - p*fvc::div(phi/fvc::interpolate(rho))
+        DpDt
     );
 
     hEqn.relax();
diff --git a/applications/solvers/heatTransfer/buoyantFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantFoam/pEqn.H
index 36378f7f6f..1cdbd0b199 100644
--- a/applications/solvers/heatTransfer/buoyantFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantFoam/pEqn.H
@@ -40,7 +40,7 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 }
 
 p == pd + rho*gh + pRef;
-dpdt = fvc::ddt(p);
+DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
 
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"

From 6c278d6fea2c3143b2271b3b0b344d4f5ed0ad38 Mon Sep 17 00:00:00 2001
From: andy <a.heather@opencfd.co.uk>
Date: Fri, 30 Jan 2009 16:18:38 +0000
Subject: [PATCH 2/5] corrected initialisation of sub-grid scale fields

---
 .../DeardorffDiffStress/DeardorffDiffStress.C | 14 ++++++--
 .../DeardorffDiffStress/DeardorffDiffStress.H |  3 ++
 .../LES/LRRDiffStress/LRRDiffStress.C         | 14 ++++++--
 .../LES/LRRDiffStress/LRRDiffStress.H         |  3 ++
 .../LES/Smagorinsky/Smagorinsky.C             | 17 ++++++---
 .../LES/Smagorinsky/Smagorinsky.H             |  3 ++
 .../LES/SpalartAllmaras/SpalartAllmaras.C     | 17 +++++++--
 .../LES/SpalartAllmaras/SpalartAllmaras.H     |  3 ++
 .../IDDESDelta/IDDESDelta.C                   |  2 +-
 .../dynMixedSmagorinsky/dynMixedSmagorinsky.C |  2 +-
 .../LES/dynOneEqEddy/dynOneEqEddy.C           | 12 +++++--
 .../LES/dynOneEqEddy/dynOneEqEddy.H           |  3 ++
 .../LES/dynSmagorinsky/dynSmagorinsky.C       | 15 +++++---
 .../LES/dynSmagorinsky/dynSmagorinsky.H       |  3 ++
 .../LES/kOmegaSSTSAS/kOmegaSSTSAS.C           | 19 +++++++---
 .../LES/kOmegaSSTSAS/kOmegaSSTSAS.H           |  6 +++-
 .../LES/locDynOneEqEddy/locDynOneEqEddy.C     | 17 +++++++--
 .../LES/locDynOneEqEddy/locDynOneEqEddy.H     |  7 ++++
 .../incompressible/LES/oneEqEddy/oneEqEddy.C  | 15 ++++++--
 .../incompressible/LES/oneEqEddy/oneEqEddy.H  |  3 ++
 .../LES/spectEddyVisc/spectEddyVisc.C         | 36 +++++++++++--------
 .../LES/spectEddyVisc/spectEddyVisc.H         |  3 ++
 22 files changed, 174 insertions(+), 43 deletions(-)

diff --git a/src/turbulenceModels/incompressible/LES/DeardorffDiffStress/DeardorffDiffStress.C b/src/turbulenceModels/incompressible/LES/DeardorffDiffStress/DeardorffDiffStress.C
index 0e3004edd7..271b7c0b92 100644
--- a/src/turbulenceModels/incompressible/LES/DeardorffDiffStress/DeardorffDiffStress.C
+++ b/src/turbulenceModels/incompressible/LES/DeardorffDiffStress/DeardorffDiffStress.C
@@ -41,6 +41,15 @@ namespace LESModels
 defineTypeNameAndDebug(DeardorffDiffStress, 0);
 addToRunTimeSelectionTable(LESModel, DeardorffDiffStress, dictionary);
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void DeardorffDiffStress::updateSubGridScaleFields(const volScalarField& K)
+{
+    nuSgs_ = ck_*sqrt(K)*delta();
+    nuSgs_.correctBoundaryConditions();
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 DeardorffDiffStress::DeardorffDiffStress
@@ -72,6 +81,8 @@ DeardorffDiffStress::DeardorffDiffStress
         )
     )
 {
+    updateSubGridScaleFields(0.5*tr(B_));
+
     printCoeffs();
 }
 
@@ -121,8 +132,7 @@ void DeardorffDiffStress::correct(const tmp<volTensorField>& tgradU)
     K = 0.5*tr(B_);
     bound(K, k0());
 
-    nuSgs_ = ck_*sqrt(K)*delta();
-    nuSgs_.correctBoundaryConditions();
+    updateSubGridScaleFields(K);
 }
 
 
diff --git a/src/turbulenceModels/incompressible/LES/DeardorffDiffStress/DeardorffDiffStress.H b/src/turbulenceModels/incompressible/LES/DeardorffDiffStress/DeardorffDiffStress.H
index bdffcf061c..67d422d64b 100644
--- a/src/turbulenceModels/incompressible/LES/DeardorffDiffStress/DeardorffDiffStress.H
+++ b/src/turbulenceModels/incompressible/LES/DeardorffDiffStress/DeardorffDiffStress.H
@@ -81,6 +81,9 @@ class DeardorffDiffStress
 
     // Private Member Functions
 
+        //- Update sub-grid scale fields
+        void updateSubGridScaleFields(const volScalarField& K);
+
         // Disallow default bitwise copy construct and assignment
         DeardorffDiffStress(const DeardorffDiffStress&);
         DeardorffDiffStress& operator=(const DeardorffDiffStress&);
diff --git a/src/turbulenceModels/incompressible/LES/LRRDiffStress/LRRDiffStress.C b/src/turbulenceModels/incompressible/LES/LRRDiffStress/LRRDiffStress.C
index 3b915025e9..3931684b6a 100644
--- a/src/turbulenceModels/incompressible/LES/LRRDiffStress/LRRDiffStress.C
+++ b/src/turbulenceModels/incompressible/LES/LRRDiffStress/LRRDiffStress.C
@@ -41,6 +41,15 @@ namespace LESModels
 defineTypeNameAndDebug(LRRDiffStress, 0);
 addToRunTimeSelectionTable(LESModel, LRRDiffStress, dictionary);
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void LRRDiffStress::updateSubGridScaleFields(const volScalarField& K)
+{
+    nuSgs_ = ck_*sqrt(K)*delta();
+    nuSgs_.correctBoundaryConditions();
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 LRRDiffStress::LRRDiffStress
@@ -81,6 +90,8 @@ LRRDiffStress::LRRDiffStress
         )
     )
 {
+    updateSubGridScaleFields(0.5*tr(B_));
+
     printCoeffs();
 }
 
@@ -131,8 +142,7 @@ void LRRDiffStress::correct(const tmp<volTensorField>& tgradU)
     K = 0.5*tr(B_);
     bound(K, k0());
 
-    nuSgs_ = ck_*sqrt(K)*delta();
-    nuSgs_.correctBoundaryConditions();
+    updateSubGridScaleFields(K);
 }
 
 
diff --git a/src/turbulenceModels/incompressible/LES/LRRDiffStress/LRRDiffStress.H b/src/turbulenceModels/incompressible/LES/LRRDiffStress/LRRDiffStress.H
index 55bb7c5d48..586d4a9f5d 100644
--- a/src/turbulenceModels/incompressible/LES/LRRDiffStress/LRRDiffStress.H
+++ b/src/turbulenceModels/incompressible/LES/LRRDiffStress/LRRDiffStress.H
@@ -80,6 +80,9 @@ class LRRDiffStress
 
     // Private Member Functions
 
+        //- Update sub-grid scale fields
+        void updateSubGridScaleFields(const volScalarField& K);
+
         // Disallow default bitwise copy construct and assignment
         LRRDiffStress(const LRRDiffStress&);
         LRRDiffStress& operator=(const LRRDiffStress&);
diff --git a/src/turbulenceModels/incompressible/LES/Smagorinsky/Smagorinsky.C b/src/turbulenceModels/incompressible/LES/Smagorinsky/Smagorinsky.C
index 0435908316..d19c2c7898 100644
--- a/src/turbulenceModels/incompressible/LES/Smagorinsky/Smagorinsky.C
+++ b/src/turbulenceModels/incompressible/LES/Smagorinsky/Smagorinsky.C
@@ -41,6 +41,15 @@ namespace LESModels
 defineTypeNameAndDebug(Smagorinsky, 0);
 addToRunTimeSelectionTable(LESModel, Smagorinsky, dictionary);
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Smagorinsky::updateSubGridScaleFields(const volTensorField& gradU)
+{
+    nuSgs_ = ck_*delta()*sqrt(k(gradU));
+    nuSgs_.correctBoundaryConditions();
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Smagorinsky::Smagorinsky
@@ -62,7 +71,9 @@ Smagorinsky::Smagorinsky
             0.094
         )
     )
-{}
+{
+    updateSubGridScaleFields(fvc::grad(U));
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -70,9 +81,7 @@ Smagorinsky::Smagorinsky
 void Smagorinsky::correct(const tmp<volTensorField>& gradU)
 {
     GenEddyVisc::correct(gradU);
-
-    nuSgs_ = ck_*delta()*sqrt(k(gradU));
-    nuSgs_.correctBoundaryConditions();
+    updateSubGridScaleFields(gradU());
 }
 
 
diff --git a/src/turbulenceModels/incompressible/LES/Smagorinsky/Smagorinsky.H b/src/turbulenceModels/incompressible/LES/Smagorinsky/Smagorinsky.H
index eba420101d..13966f4883 100644
--- a/src/turbulenceModels/incompressible/LES/Smagorinsky/Smagorinsky.H
+++ b/src/turbulenceModels/incompressible/LES/Smagorinsky/Smagorinsky.H
@@ -77,6 +77,9 @@ class Smagorinsky
 
     // Private Member Functions
 
+        //- Update sub-grid scale fields
+        void updateSubGridScaleFields(const volTensorField& gradU);
+
         // Disallow default bitwise copy construct and assignment
         Smagorinsky(const Smagorinsky&);
         Smagorinsky& operator=(const Smagorinsky&);
diff --git a/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.C b/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.C
index a6d7b44bb5..91ee6900b0 100644
--- a/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.C
+++ b/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.C
@@ -41,6 +41,15 @@ namespace LESModels
 defineTypeNameAndDebug(SpalartAllmaras, 0);
 addToRunTimeSelectionTable(LESModel, SpalartAllmaras, dictionary);
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void SpalartAllmaras::updateSubGridScaleFields()
+{
+    nuSgs_.internalField() = fv1()*nuTilda_.internalField();
+    nuSgs_.correctBoundaryConditions();
+}
+
+
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 tmp<volScalarField> SpalartAllmaras::fv1() const
@@ -266,7 +275,10 @@ SpalartAllmaras::SpalartAllmaras
         ),
         mesh_
     )
-{}
+{
+    updateSubGridScaleFields();
+}
+
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -306,8 +318,7 @@ void SpalartAllmaras::correct(const tmp<volTensorField>& gradU)
     bound(nuTilda_, dimensionedScalar("zero", nuTilda_.dimensions(), 0.0));
     nuTilda_.correctBoundaryConditions();
 
-    nuSgs_.internalField() = fv1()*nuTilda_.internalField();
-    nuSgs_.correctBoundaryConditions();
+    updateSubGridScaleFields();
 }
 
 
diff --git a/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.H b/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.H
index b813075b1d..60dfbcc555 100644
--- a/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.H
+++ b/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.H
@@ -59,6 +59,9 @@ class SpalartAllmaras
 {
     // Private member functions
 
+        //- Update sub-grid scale fields
+        void updateSubGridScaleFields();
+
         // Disallow default bitwise copy construct and assignment
         SpalartAllmaras(const SpalartAllmaras&);
         SpalartAllmaras& operator=(const SpalartAllmaras&);
diff --git a/src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/IDDESDelta/IDDESDelta.C b/src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/IDDESDelta/IDDESDelta.C
index ef02e48456..3f3aa7699b 100644
--- a/src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/IDDESDelta/IDDESDelta.C
+++ b/src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/IDDESDelta/IDDESDelta.C
@@ -44,7 +44,7 @@ void Foam::IDDESDelta::calcDelta()
     const Vector<label>& directions = mesh().directions();
     label nD = (directions.nComponents + cmptSum(directions))/2;
 
-    // - Init hwn as wall distant.
+    // initialise hwn as wall distance
     volScalarField hwn = wallDist(mesh()).y();
 
     scalar deltamaxTmp = 0.;
diff --git a/src/turbulenceModels/incompressible/LES/dynMixedSmagorinsky/dynMixedSmagorinsky.C b/src/turbulenceModels/incompressible/LES/dynMixedSmagorinsky/dynMixedSmagorinsky.C
index c0fb5160d3..43f684d228 100644
--- a/src/turbulenceModels/incompressible/LES/dynMixedSmagorinsky/dynMixedSmagorinsky.C
+++ b/src/turbulenceModels/incompressible/LES/dynMixedSmagorinsky/dynMixedSmagorinsky.C
@@ -113,7 +113,7 @@ tmp<fvVectorMatrix> dynMixedSmagorinsky::divDevBeff(volVectorField& U) const
 void dynMixedSmagorinsky::correct(const tmp<volTensorField>& gradU)
 {
     scaleSimilarity::correct(gradU);
-    dynSmagorinsky::correct(gradU);
+    dynSmagorinsky::correct(gradU());
 }
 
 
diff --git a/src/turbulenceModels/incompressible/LES/dynOneEqEddy/dynOneEqEddy.C b/src/turbulenceModels/incompressible/LES/dynOneEqEddy/dynOneEqEddy.C
index e447211bec..4342712d80 100644
--- a/src/turbulenceModels/incompressible/LES/dynOneEqEddy/dynOneEqEddy.C
+++ b/src/turbulenceModels/incompressible/LES/dynOneEqEddy/dynOneEqEddy.C
@@ -43,6 +43,13 @@ addToRunTimeSelectionTable(LESModel, dynOneEqEddy, dictionary);
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+void dynOneEqEddy::updateSubGridScaleFields(const volSymmTensorField& D)
+{
+    nuSgs_ = ck(D)*sqrt(k_)*delta();
+    nuSgs_.correctBoundaryConditions();
+}
+
+
 dimensionedScalar dynOneEqEddy::ck(const volSymmTensorField& D) const
 {
     volScalarField KK = 0.5*(filter_(magSqr(U())) - magSqr(filter_(U())));
@@ -120,6 +127,8 @@ dynOneEqEddy::dynOneEqEddy
     filterPtr_(LESfilter::New(U.mesh(), coeffDict())),
     filter_(filterPtr_())
 {
+    updateSubGridScaleFields(symm(fvc::grad(U)));
+
     printCoeffs();
 }
 
@@ -155,8 +164,7 @@ void dynOneEqEddy::correct(const tmp<volTensorField>& gradU)
 
     bound(k_, k0());
 
-    nuSgs_ = ck(D)*sqrt(k_)*delta();
-    nuSgs_.correctBoundaryConditions();
+    updateSubGridScaleFields(D);
 }
 
 
diff --git a/src/turbulenceModels/incompressible/LES/dynOneEqEddy/dynOneEqEddy.H b/src/turbulenceModels/incompressible/LES/dynOneEqEddy/dynOneEqEddy.H
index b8cb1da312..17ed2c62ad 100644
--- a/src/turbulenceModels/incompressible/LES/dynOneEqEddy/dynOneEqEddy.H
+++ b/src/turbulenceModels/incompressible/LES/dynOneEqEddy/dynOneEqEddy.H
@@ -87,6 +87,9 @@ class dynOneEqEddy
 
     // Private Member Functions
 
+        //- Update sub-grid scale fields
+        void updateSubGridScaleFields(const volSymmTensorField& D);
+
         //- Calculate ck, ce by filtering the velocity field U.
         dimensionedScalar ck(const volSymmTensorField& D) const;
         dimensionedScalar ce(const volSymmTensorField& D) const;
diff --git a/src/turbulenceModels/incompressible/LES/dynSmagorinsky/dynSmagorinsky.C b/src/turbulenceModels/incompressible/LES/dynSmagorinsky/dynSmagorinsky.C
index 6d1f711956..e260a6d16d 100644
--- a/src/turbulenceModels/incompressible/LES/dynSmagorinsky/dynSmagorinsky.C
+++ b/src/turbulenceModels/incompressible/LES/dynSmagorinsky/dynSmagorinsky.C
@@ -43,6 +43,13 @@ addToRunTimeSelectionTable(LESModel, dynSmagorinsky, dictionary);
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+void dynSmagorinsky::updateSubGridScaleFields(const volSymmTensorField& D)
+{
+    nuSgs_ = cD(D)*sqr(delta())*sqrt(magSqr(D));
+    nuSgs_.correctBoundaryConditions();
+}
+
+
 dimensionedScalar dynSmagorinsky::cD(const volSymmTensorField& D) const
 {
     volSymmTensorField LL = dev(filter_(sqr(U())) - (sqr(filter_(U()))));
@@ -111,6 +118,8 @@ dynSmagorinsky::dynSmagorinsky
     filterPtr_(LESfilter::New(U.mesh(), coeffDict())),
     filter_(filterPtr_())
 {
+    updateSubGridScaleFields(dev(symm(fvc::grad(U))));
+
     printCoeffs();
 }
 
@@ -128,12 +137,10 @@ void dynSmagorinsky::correct(const tmp<volTensorField>& gradU)
     LESModel::correct(gradU);
 
     volSymmTensorField D = dev(symm(gradU));
-    volScalarField magSqrD = magSqr(D);
 
-    k_ = cI(D)*sqr(delta())*magSqrD;
+    k_ = cI(D)*sqr(delta())*magSqr(D);
 
-    nuSgs_ = cD(D)*sqr(delta())*sqrt(magSqrD);
-    nuSgs_.correctBoundaryConditions();
+    updateSubGridScaleFields(D);
 }
 
 
diff --git a/src/turbulenceModels/incompressible/LES/dynSmagorinsky/dynSmagorinsky.H b/src/turbulenceModels/incompressible/LES/dynSmagorinsky/dynSmagorinsky.H
index 2fa29ff41a..39a6636572 100644
--- a/src/turbulenceModels/incompressible/LES/dynSmagorinsky/dynSmagorinsky.H
+++ b/src/turbulenceModels/incompressible/LES/dynSmagorinsky/dynSmagorinsky.H
@@ -96,6 +96,9 @@ class dynSmagorinsky
 
     // Private Member Functions
 
+        //- Update sub-grid scale fields
+        void updateSubGridScaleFields(const volSymmTensorField& D);
+
         //- Calculate coefficients cD, cI from filtering velocity field
         dimensionedScalar cD(const volSymmTensorField& D) const;
         dimensionedScalar cI(const volSymmTensorField& D) const;
diff --git a/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C b/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C
index 3460137709..b4aca6c7f3 100644
--- a/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C
+++ b/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C
@@ -44,6 +44,13 @@ addToRunTimeSelectionTable(LESModel, kOmegaSSTSAS, dictionary);
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
+void kOmegaSSTSAS::updateSubGridScaleFields(const volScalarField& S2)
+{
+    nuSgs_ == a1_*k_/max(a1_*omega_, F2()*sqrt(S2));
+    nuSgs_.correctBoundaryConditions();
+}
+
+
 tmp<volScalarField> kOmegaSSTSAS::F1(const volScalarField& CDkOmega) const
 {
     volScalarField CDkOmegaPlus = max
@@ -304,6 +311,8 @@ kOmegaSSTSAS::kOmegaSSTSAS
         mesh_
     )
 {
+    updateSubGridScaleFields(magSqr(symm(fvc::grad(U))));
+
     printCoeffs();
 }
 
@@ -325,7 +334,8 @@ void kOmegaSSTSAS::correct(const tmp<volTensorField>& gradU)
     volVectorField gradK = fvc::grad(k_);
     volVectorField gradOmega = fvc::grad(omega_);
     volScalarField L = sqrt(k_)/(pow(Cmu_, 0.25)*(omega_ + omegaSmall_));
-    volScalarField CDkOmega = (2.0*alphaOmega2_)*(gradK & gradOmega)/(omega_ + omegaSmall_);
+    volScalarField CDkOmega =
+        (2.0*alphaOmega2_)*(gradK & gradOmega)/(omega_ + omegaSmall_);
     volScalarField F1 = this->F1(CDkOmega);
     volScalarField G = nuSgs_*2.0*S2;
 
@@ -375,7 +385,8 @@ void kOmegaSSTSAS::correct(const tmp<volTensorField>& gradU)
            *max
             (
                 dimensionedScalar("zero",dimensionSet(0, 0 , -2, 0, 0),0. ),
-                zetaTilda2_*kappa_*S2*(L/Lvk2(S2))- 2.0/alphaPhi_*k_*grad_omega_k
+                zetaTilda2_*kappa_*S2*(L/Lvk2(S2))
+              - 2.0/alphaPhi_*k_*grad_omega_k
             )
         );
 
@@ -384,9 +395,7 @@ void kOmegaSSTSAS::correct(const tmp<volTensorField>& gradU)
     }
     bound(omega_, omega0_);
 
-    // Re-calculate viscosity
-    nuSgs_ == a1_*k_/max(a1_*omega_, F2()*sqrt(S2));
-    nuSgs_.correctBoundaryConditions();
+    updateSubGridScaleFields(S2);
 }
 
 
diff --git a/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.H b/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.H
index 56e395ff4c..e5b03aef68 100644
--- a/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.H
+++ b/src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.H
@@ -62,6 +62,9 @@ class kOmegaSSTSAS
 {
     // Private member functions
 
+        //- Update sub-grid scale fields
+        void updateSubGridScaleFields(const volScalarField& D);
+
         // Disallow default bitwise copy construct and assignment
         kOmegaSSTSAS(const kOmegaSSTSAS&);
         kOmegaSSTSAS& operator=(const kOmegaSSTSAS&);
@@ -237,7 +240,8 @@ public:
         //  i.e. the additional term in the filtered NSE.
         virtual tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const;
 
-        //- Solve the turbulence equations (k-w) and correct the turbulence viscosity
+        //- Solve the turbulence equations (k-w) and correct the turbulence
+        //  viscosity
         virtual void correct(const tmp<volTensorField>& gradU);
 
         //- Read turbulenceProperties dictionary
diff --git a/src/turbulenceModels/incompressible/LES/locDynOneEqEddy/locDynOneEqEddy.C b/src/turbulenceModels/incompressible/LES/locDynOneEqEddy/locDynOneEqEddy.C
index a195723af7..7d34bbb981 100644
--- a/src/turbulenceModels/incompressible/LES/locDynOneEqEddy/locDynOneEqEddy.C
+++ b/src/turbulenceModels/incompressible/LES/locDynOneEqEddy/locDynOneEqEddy.C
@@ -43,6 +43,17 @@ addToRunTimeSelectionTable(LESModel, locDynOneEqEddy, dictionary);
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+void locDynOneEqEddy::updateSubGridScaleFields
+(
+    const volSymmTensorField& D,
+    const volScalarField& KK
+)
+{
+    nuSgs_ = ck(D, KK)*sqrt(k_)*delta();
+    nuSgs_.correctBoundaryConditions();
+}
+
+
 volScalarField locDynOneEqEddy::ck
 (
     const volSymmTensorField& D,
@@ -108,6 +119,9 @@ locDynOneEqEddy::locDynOneEqEddy
     filterPtr_(LESfilter::New(U.mesh(), coeffDict())),
     filter_(filterPtr_())
 {
+    volScalarField KK = 0.5*(filter_(magSqr(U)) - magSqr(filter_(U)));
+    updateSubGridScaleFields(symm(fvc::grad(U)), KK);
+
     printCoeffs();
 }
 
@@ -146,8 +160,7 @@ void locDynOneEqEddy::correct(const tmp<volTensorField>& gradU)
 
     bound(k_, k0());
 
-    nuSgs_ = ck(D, KK)*sqrt(k_)*delta();
-    nuSgs_.correctBoundaryConditions();
+    updateSubGridScaleFields(D, KK);
 }
 
 
diff --git a/src/turbulenceModels/incompressible/LES/locDynOneEqEddy/locDynOneEqEddy.H b/src/turbulenceModels/incompressible/LES/locDynOneEqEddy/locDynOneEqEddy.H
index 6b142c4e37..4848c6492e 100644
--- a/src/turbulenceModels/incompressible/LES/locDynOneEqEddy/locDynOneEqEddy.H
+++ b/src/turbulenceModels/incompressible/LES/locDynOneEqEddy/locDynOneEqEddy.H
@@ -96,6 +96,13 @@ class locDynOneEqEddy
 
     // Private Member Functions
 
+        //- Update sub-grid scale fields
+        void updateSubGridScaleFields
+        (
+            const volSymmTensorField& D,
+            const volScalarField& KK
+        );
+
         //- Calculate ck, ce by filtering the velocity field U.
         volScalarField ck
         (
diff --git a/src/turbulenceModels/incompressible/LES/oneEqEddy/oneEqEddy.C b/src/turbulenceModels/incompressible/LES/oneEqEddy/oneEqEddy.C
index da55eeb949..d6a5de4ed4 100644
--- a/src/turbulenceModels/incompressible/LES/oneEqEddy/oneEqEddy.C
+++ b/src/turbulenceModels/incompressible/LES/oneEqEddy/oneEqEddy.C
@@ -41,6 +41,16 @@ namespace LESModels
 defineTypeNameAndDebug(oneEqEddy, 0);
 addToRunTimeSelectionTable(LESModel, oneEqEddy, dictionary);
 
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void oneEqEddy::updateSubGridScaleFields()
+{
+    nuSgs_ = ck_*sqrt(k_)*delta();
+    nuSgs_.correctBoundaryConditions();
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 oneEqEddy::oneEqEddy
@@ -76,6 +86,8 @@ oneEqEddy::oneEqEddy
         )
     )
 {
+    updateSubGridScaleFields();
+
     printCoeffs();
 }
 
@@ -103,8 +115,7 @@ void oneEqEddy::correct(const tmp<volTensorField>& gradU)
 
     bound(k_, k0());
 
-    nuSgs_ = ck_*sqrt(k_)*delta();
-    nuSgs_.correctBoundaryConditions();
+    updateSubGridScaleFields();
 }
 
 
diff --git a/src/turbulenceModels/incompressible/LES/oneEqEddy/oneEqEddy.H b/src/turbulenceModels/incompressible/LES/oneEqEddy/oneEqEddy.H
index 65b01ca479..c060f82d67 100644
--- a/src/turbulenceModels/incompressible/LES/oneEqEddy/oneEqEddy.H
+++ b/src/turbulenceModels/incompressible/LES/oneEqEddy/oneEqEddy.H
@@ -83,6 +83,9 @@ class oneEqEddy
 
     // Private Member Functions
 
+        //- Update sub-grid scale fields
+        void updateSubGridScaleFields();
+
         // Disallow default bitwise copy construct and assignment
         oneEqEddy(const oneEqEddy&);
         oneEqEddy& operator=(const oneEqEddy&);
diff --git a/src/turbulenceModels/incompressible/LES/spectEddyVisc/spectEddyVisc.C b/src/turbulenceModels/incompressible/LES/spectEddyVisc/spectEddyVisc.C
index 28a563b4c6..4d8a7edc44 100644
--- a/src/turbulenceModels/incompressible/LES/spectEddyVisc/spectEddyVisc.C
+++ b/src/turbulenceModels/incompressible/LES/spectEddyVisc/spectEddyVisc.C
@@ -41,6 +41,25 @@ namespace LESModels
 defineTypeNameAndDebug(spectEddyVisc, 0);
 addToRunTimeSelectionTable(LESModel, spectEddyVisc, dictionary);
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void spectEddyVisc::updateSubGridScaleFields(const volTensorField& gradU)
+{
+    volScalarField Re = sqr(delta())*mag(symm(gradU))/nu();
+    for (label i=0; i<5; i++)
+    {
+        nuSgs_ =
+            nu()
+           /(
+                 scalar(1)
+               - exp(-cB_*pow(nu()/(nuSgs_ + nu()), 1.0/3.0)*pow(Re, -2.0/3.0))
+            );
+    }
+
+    nuSgs_.correctBoundaryConditions();
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 spectEddyVisc::spectEddyVisc
@@ -101,6 +120,8 @@ spectEddyVisc::spectEddyVisc
     )
 {
     printCoeffs();
+
+    updateSubGridScaleFields(fvc::grad(U));
 }
 
 
@@ -121,20 +142,7 @@ tmp<volScalarField> spectEddyVisc::k() const
 void spectEddyVisc::correct(const tmp<volTensorField>& gradU)
 {
     GenEddyVisc::correct(gradU);
-
-    volScalarField Re = sqr(delta())*mag(symm(gradU))/nu();
-
-    for (label i=0; i<5; i++)
-    {
-        nuSgs_ =
-            nu()
-           /(
-               scalar(1)
-               - exp(-cB_*pow(nu()/(nuSgs_ + nu()), 1.0/3.0)*pow(Re, -2.0/3.0))
-            );
-    }
-
-    nuSgs_.correctBoundaryConditions();
+    updateSubGridScaleFields(gradU());
 }
 
 
diff --git a/src/turbulenceModels/incompressible/LES/spectEddyVisc/spectEddyVisc.H b/src/turbulenceModels/incompressible/LES/spectEddyVisc/spectEddyVisc.H
index 715b65c4b5..709c9ef07c 100644
--- a/src/turbulenceModels/incompressible/LES/spectEddyVisc/spectEddyVisc.H
+++ b/src/turbulenceModels/incompressible/LES/spectEddyVisc/spectEddyVisc.H
@@ -85,6 +85,9 @@ class spectEddyVisc
 
     // Private Member Functions
 
+        //- Update sub-grid scale fields
+        void updateSubGridScaleFields(const volTensorField& gradU);
+
         // Disallow default bitwise copy construct and assignment
         spectEddyVisc(const spectEddyVisc&);
         spectEddyVisc& operator=(const spectEddyVisc&);

From 534d06054e15c1a5ca49552213db5bb10e074722 Mon Sep 17 00:00:00 2001
From: andy <a.heather@opencfd.co.uk>
Date: Fri, 30 Jan 2009 16:37:03 +0000
Subject: [PATCH 3/5] corrected initialisation of sub-grid scale fields -
 compressible

---
 .../DeardorffDiffStress/DeardorffDiffStress.C | 15 +++++++--
 .../DeardorffDiffStress/DeardorffDiffStress.H |  3 ++
 .../LES/Smagorinsky/Smagorinsky.C             | 32 ++++++++++++-------
 .../LES/Smagorinsky/Smagorinsky.H             |  3 ++
 .../LES/SpalartAllmaras/SpalartAllmaras.C     | 14 ++++++--
 .../LES/SpalartAllmaras/SpalartAllmaras.H     |  5 ++-
 .../LES/dynOneEqEddy/dynOneEqEddy.C           | 12 +++++--
 .../LES/dynOneEqEddy/dynOneEqEddy.H           |  3 ++
 .../LES/lowReOneEqEddy/lowReOneEqEddy.C       | 23 +++++++++----
 .../LES/lowReOneEqEddy/lowReOneEqEddy.H       |  3 ++
 .../compressible/LES/oneEqEddy/oneEqEddy.C    | 14 ++++++--
 .../compressible/LES/oneEqEddy/oneEqEddy.H    |  3 ++
 12 files changed, 102 insertions(+), 28 deletions(-)

diff --git a/src/turbulenceModels/compressible/LES/DeardorffDiffStress/DeardorffDiffStress.C b/src/turbulenceModels/compressible/LES/DeardorffDiffStress/DeardorffDiffStress.C
index aa0b474992..e847acda6e 100644
--- a/src/turbulenceModels/compressible/LES/DeardorffDiffStress/DeardorffDiffStress.C
+++ b/src/turbulenceModels/compressible/LES/DeardorffDiffStress/DeardorffDiffStress.C
@@ -41,6 +41,16 @@ namespace LESModels
 defineTypeNameAndDebug(DeardorffDiffStress, 0);
 addToRunTimeSelectionTable(LESModel, DeardorffDiffStress, dictionary);
 
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void DeardorffDiffStress::updateSubGridScaleFields(const volScalarField& K)
+{
+    muSgs_ = ck_*rho()*sqrt(K)*delta();
+    muSgs_.correctBoundaryConditions();
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 DeardorffDiffStress::DeardorffDiffStress
@@ -73,6 +83,8 @@ DeardorffDiffStress::DeardorffDiffStress
         )
     )
 {
+    updateSubGridScaleFields(0.5*tr(B_));
+
     printCoeffs();
 }
 
@@ -119,8 +131,7 @@ void DeardorffDiffStress::correct(const tmp<volTensorField>& tgradU)
     K = 0.5*tr(B_);
     bound(K, k0());
 
-    muSgs_ = ck_*rho()*sqrt(K)*delta();
-    muSgs_.correctBoundaryConditions();
+    updateSubGridScaleFields(K);
 }
 
 
diff --git a/src/turbulenceModels/compressible/LES/DeardorffDiffStress/DeardorffDiffStress.H b/src/turbulenceModels/compressible/LES/DeardorffDiffStress/DeardorffDiffStress.H
index 4371afbc66..510c304f20 100644
--- a/src/turbulenceModels/compressible/LES/DeardorffDiffStress/DeardorffDiffStress.H
+++ b/src/turbulenceModels/compressible/LES/DeardorffDiffStress/DeardorffDiffStress.H
@@ -81,6 +81,9 @@ class DeardorffDiffStress
 
     // Private Member Functions
 
+        //- Update sub-grid scale fields
+        void updateSubGridScaleFields(const volScalarField& K);
+
         // Disallow default bitwise copy construct and assignment
         DeardorffDiffStress(const DeardorffDiffStress&);
         DeardorffDiffStress& operator=(const DeardorffDiffStress&);
diff --git a/src/turbulenceModels/compressible/LES/Smagorinsky/Smagorinsky.C b/src/turbulenceModels/compressible/LES/Smagorinsky/Smagorinsky.C
index e845776031..c89b2aeecf 100644
--- a/src/turbulenceModels/compressible/LES/Smagorinsky/Smagorinsky.C
+++ b/src/turbulenceModels/compressible/LES/Smagorinsky/Smagorinsky.C
@@ -41,6 +41,24 @@ namespace LESModels
 defineTypeNameAndDebug(Smagorinsky, 0);
 addToRunTimeSelectionTable(LESModel, Smagorinsky, dictionary);
 
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Smagorinsky::updateSubGridScaleFields(const volTensorField& gradU)
+{
+    volSymmTensorField D = symm(gradU);
+
+    volScalarField a = ce_/delta();
+    volScalarField b = (2.0/3.0)*tr(D);
+    volScalarField c = 2*ck_*delta()*(dev(D) && D);
+
+    k_ = sqr((-b + sqrt(sqr(b) + 4*a*c))/(2*a));
+
+    muSgs_ = ck_*rho()*delta()*sqrt(k_);
+    muSgs_.correctBoundaryConditions();
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Smagorinsky::Smagorinsky
@@ -64,6 +82,8 @@ Smagorinsky::Smagorinsky
         )
     )
 {
+    updateSubGridScaleFields(fvc::grad(U));
+
     printCoeffs();
 }
 
@@ -73,17 +93,7 @@ Smagorinsky::Smagorinsky
 void Smagorinsky::correct(const tmp<volTensorField>& gradU)
 {
     GenEddyVisc::correct(gradU);
-
-    volSymmTensorField D = symm(gradU);
-
-    volScalarField a = ce_/delta();
-    volScalarField b = (2.0/3.0)*tr(D);
-    volScalarField c = 2*ck_*delta()*(dev(D) && D);
-
-    k_ = sqr((-b + sqrt(sqr(b) + 4*a*c))/(2*a));
-
-    muSgs_ = ck_*rho()*delta()*sqrt(k_);
-    muSgs_.correctBoundaryConditions();
+    updateSubGridScaleFields(gradU());
 }
 
 
diff --git a/src/turbulenceModels/compressible/LES/Smagorinsky/Smagorinsky.H b/src/turbulenceModels/compressible/LES/Smagorinsky/Smagorinsky.H
index c93745dbcb..ba296a87f9 100644
--- a/src/turbulenceModels/compressible/LES/Smagorinsky/Smagorinsky.H
+++ b/src/turbulenceModels/compressible/LES/Smagorinsky/Smagorinsky.H
@@ -76,6 +76,9 @@ class Smagorinsky
 
     // Private Member Functions
 
+        //- Update sub-grid scale fields
+        void updateSubGridScaleFields(const volTensorField& gradU);
+
         // Disallow default bitwise copy construct and assignment
         Smagorinsky(const Smagorinsky&);
         Smagorinsky& operator=(const Smagorinsky&);
diff --git a/src/turbulenceModels/compressible/LES/SpalartAllmaras/SpalartAllmaras.C b/src/turbulenceModels/compressible/LES/SpalartAllmaras/SpalartAllmaras.C
index 92d15f5cd8..b8d1a8b97a 100644
--- a/src/turbulenceModels/compressible/LES/SpalartAllmaras/SpalartAllmaras.C
+++ b/src/turbulenceModels/compressible/LES/SpalartAllmaras/SpalartAllmaras.C
@@ -45,6 +45,13 @@ addToRunTimeSelectionTable(LESModel, SpalartAllmaras, dictionary);
 
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
+void SpalartAllmaras::updateSubGridScaleFields()
+{
+    muSgs_.internalField() = rho()*fv1()*nuTilda_.internalField();
+    muSgs_.correctBoundaryConditions();
+}
+
+
 tmp<volScalarField> SpalartAllmaras::fv1() const
 {
     volScalarField chi3 = pow3(nuTilda_/(mu()/rho()));
@@ -225,6 +232,8 @@ SpalartAllmaras::SpalartAllmaras
     )
 
 {
+    updateSubGridScaleFields();
+
     printCoeffs();
 }
 
@@ -288,10 +297,9 @@ void SpalartAllmaras::correct(const tmp<volTensorField>& tgradU)
     );
 
     bound(nuTilda_, dimensionedScalar("zero", nuTilda_.dimensions(), 0.0));
-
     nuTilda_.correctBoundaryConditions();
-    muSgs_.internalField() = rho()*fv1()*nuTilda_.internalField();
-    muSgs_.correctBoundaryConditions();
+
+    updateSubGridScaleFields();
 }
 
 
diff --git a/src/turbulenceModels/compressible/LES/SpalartAllmaras/SpalartAllmaras.H b/src/turbulenceModels/compressible/LES/SpalartAllmaras/SpalartAllmaras.H
index 437baedb49..e119768537 100644
--- a/src/turbulenceModels/compressible/LES/SpalartAllmaras/SpalartAllmaras.H
+++ b/src/turbulenceModels/compressible/LES/SpalartAllmaras/SpalartAllmaras.H
@@ -49,7 +49,7 @@ namespace LESModels
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class SpalartAllmaras Declaration
+                       Class SpalartAllmaras Declaration
 \*---------------------------------------------------------------------------*/
 
 class SpalartAllmaras
@@ -74,6 +74,9 @@ class SpalartAllmaras
 
     // Private member functions
 
+        //- Update sub-grid scale fields
+        void updateSubGridScaleFields();
+
         tmp<volScalarField> fv1() const;
         tmp<volScalarField> fv2() const;
         tmp<volScalarField> fv3() const;
diff --git a/src/turbulenceModels/compressible/LES/dynOneEqEddy/dynOneEqEddy.C b/src/turbulenceModels/compressible/LES/dynOneEqEddy/dynOneEqEddy.C
index f6194e4882..7e9b3d9b7b 100644
--- a/src/turbulenceModels/compressible/LES/dynOneEqEddy/dynOneEqEddy.C
+++ b/src/turbulenceModels/compressible/LES/dynOneEqEddy/dynOneEqEddy.C
@@ -43,6 +43,13 @@ addToRunTimeSelectionTable(LESModel, dynOneEqEddy, dictionary);
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+void dynOneEqEddy::updateSubGridScaleFields(const volSymmTensorField& D)
+{
+    muSgs_ = ck_(D)*rho()*sqrt(k_)*delta();
+    muSgs_.correctBoundaryConditions();
+}
+
+
 dimensionedScalar dynOneEqEddy::ck_(const volSymmTensorField& D) const
 {
     volScalarField KK = 0.5*(filter_(magSqr(U())) - magSqr(filter_(U())));
@@ -90,6 +97,8 @@ dynOneEqEddy::dynOneEqEddy
     filterPtr_(LESfilter::New(U.mesh(), coeffDict())),
     filter_(filterPtr_())
 {
+    updateSubGridScaleFields(dev(symm(fvc::grad(U))));
+
     printCoeffs();
 }
 
@@ -119,8 +128,7 @@ void dynOneEqEddy::correct(const tmp<volTensorField>& tgradU)
 
     bound(k_, dimensionedScalar("0", k_.dimensions(), 1.0e-10));
 
-    muSgs_ = ck_(D)*rho()*sqrt(k_)*delta();
-    muSgs_.correctBoundaryConditions();
+    updateSubGridScaleFields(D);
 }
 
 
diff --git a/src/turbulenceModels/compressible/LES/dynOneEqEddy/dynOneEqEddy.H b/src/turbulenceModels/compressible/LES/dynOneEqEddy/dynOneEqEddy.H
index 4f5ec19b70..a32dacddfd 100644
--- a/src/turbulenceModels/compressible/LES/dynOneEqEddy/dynOneEqEddy.H
+++ b/src/turbulenceModels/compressible/LES/dynOneEqEddy/dynOneEqEddy.H
@@ -82,6 +82,9 @@ class dynOneEqEddy
 
     // Private Member Functions
 
+        //- Update sub-grid scale fields
+        void updateSubGridScaleFields(const volSymmTensorField& D);
+
         //- Calculate ck, ce by filtering the velocity field U.
         dimensionedScalar ck_(const volSymmTensorField& D) const;
         dimensionedScalar ce_(const volSymmTensorField& D) const;
diff --git a/src/turbulenceModels/compressible/LES/lowReOneEqEddy/lowReOneEqEddy.C b/src/turbulenceModels/compressible/LES/lowReOneEqEddy/lowReOneEqEddy.C
index 3783677eaf..b5670ad245 100644
--- a/src/turbulenceModels/compressible/LES/lowReOneEqEddy/lowReOneEqEddy.C
+++ b/src/turbulenceModels/compressible/LES/lowReOneEqEddy/lowReOneEqEddy.C
@@ -41,6 +41,19 @@ namespace LESModels
 defineTypeNameAndDebug(lowReOneEqEddy, 0);
 addToRunTimeSelectionTable(LESModel, lowReOneEqEddy, dictionary);
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void lowReOneEqEddy::updateSubGridScaleFields()
+{
+    // High Re eddy viscosity
+    muSgs_ = ck_*rho()*sqrt(k_)*delta();
+
+    // low Re no corrected eddy viscosity
+    muSgs_ -= (mu()/beta_)*(scalar(1) - exp(-beta_*muSgs_/mu()));
+    muSgs_.correctBoundaryConditions();
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 lowReOneEqEddy::lowReOneEqEddy
@@ -73,6 +86,8 @@ lowReOneEqEddy::lowReOneEqEddy
         )
     )
 {
+    updateSubGridScaleFields();
+
     printCoeffs();
 }
 
@@ -101,13 +116,7 @@ void lowReOneEqEddy::correct(const tmp<volTensorField>& tgradU)
 
     bound(k_, k0());
 
-    // High Re eddy viscosity
-    muSgs_ = ck_*rho()*sqrt(k_)*delta();
-
-    // low Re no corrected eddy viscosity
-    muSgs_ -= (mu()/beta_)*(scalar(1) - exp(-beta_*muSgs_/mu()));
-
-    muSgs_.correctBoundaryConditions();
+    updateSubGridScaleFields();
 }
 
 
diff --git a/src/turbulenceModels/compressible/LES/lowReOneEqEddy/lowReOneEqEddy.H b/src/turbulenceModels/compressible/LES/lowReOneEqEddy/lowReOneEqEddy.H
index a36ec44860..9820905a41 100644
--- a/src/turbulenceModels/compressible/LES/lowReOneEqEddy/lowReOneEqEddy.H
+++ b/src/turbulenceModels/compressible/LES/lowReOneEqEddy/lowReOneEqEddy.H
@@ -77,6 +77,9 @@ class lowReOneEqEddy
 
     // Private Member Functions
 
+        //- Update sub-grid scale fields
+        void updateSubGridScaleFields();
+
         // Disallow default bitwise copy construct and assignment
         lowReOneEqEddy(const lowReOneEqEddy&);
         lowReOneEqEddy& operator=(const lowReOneEqEddy&);
diff --git a/src/turbulenceModels/compressible/LES/oneEqEddy/oneEqEddy.C b/src/turbulenceModels/compressible/LES/oneEqEddy/oneEqEddy.C
index 65457c0359..414aba5b03 100644
--- a/src/turbulenceModels/compressible/LES/oneEqEddy/oneEqEddy.C
+++ b/src/turbulenceModels/compressible/LES/oneEqEddy/oneEqEddy.C
@@ -41,6 +41,15 @@ namespace LESModels
 defineTypeNameAndDebug(oneEqEddy, 0);
 addToRunTimeSelectionTable(LESModel, oneEqEddy, dictionary);
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void oneEqEddy::updateSubGridScaleFields()
+{
+    muSgs_ = ck_*rho()*sqrt(k_)*delta();
+    muSgs_.correctBoundaryConditions();
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 oneEqEddy::oneEqEddy
@@ -64,6 +73,8 @@ oneEqEddy::oneEqEddy
         )
     )
 {
+    updateSubGridScaleFields();
+
     printCoeffs();
 }
 
@@ -95,8 +106,7 @@ void oneEqEddy::correct(const tmp<volTensorField>& tgradU)
 
     bound(k_, k0());
 
-    muSgs_ = ck_*rho()*sqrt(k_)*delta();
-    muSgs_.correctBoundaryConditions();
+    updateSubGridScaleFields();
 }
 
 
diff --git a/src/turbulenceModels/compressible/LES/oneEqEddy/oneEqEddy.H b/src/turbulenceModels/compressible/LES/oneEqEddy/oneEqEddy.H
index 6568724fd2..2a9bb74591 100644
--- a/src/turbulenceModels/compressible/LES/oneEqEddy/oneEqEddy.H
+++ b/src/turbulenceModels/compressible/LES/oneEqEddy/oneEqEddy.H
@@ -80,6 +80,9 @@ class oneEqEddy
 
     // Private Member Functions
 
+        //- Update sub-grid scale fields
+        void updateSubGridScaleFields();
+
         // Disallow default bitwise copy construct and assignment
         oneEqEddy(const oneEqEddy&);
         oneEqEddy& operator=(const oneEqEddy&);

From 5afc0db5a82c384821f3e51fbe997ef052490404 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Mon, 2 Feb 2009 12:23:41 +0000
Subject: [PATCH 4/5] Solaris port

---
 etc/bashrc                          | 11 +++++++++++
 etc/cshrc                           | 13 ++++++++++++-
 etc/settings.csh                    |  9 +++++++++
 etc/settings.sh                     | 10 ++++++++++
 src/OSspecific/Unix/Make/files      |  6 +++++-
 wmake/rules/SunOS64Gcc/X            |  3 +++
 wmake/rules/SunOS64Gcc/c            | 16 ++++++++++++++++
 wmake/rules/SunOS64Gcc/c++          | 21 +++++++++++++++++++++
 wmake/rules/SunOS64Gcc/c++Debug     |  2 ++
 wmake/rules/SunOS64Gcc/c++Opt       |  2 ++
 wmake/rules/SunOS64Gcc/c++Prof      |  2 ++
 wmake/rules/SunOS64Gcc/cDebug       |  2 ++
 wmake/rules/SunOS64Gcc/cOpt         |  2 ++
 wmake/rules/SunOS64Gcc/cProf        |  2 ++
 wmake/rules/SunOS64Gcc/general      | 11 +++++++++++
 wmake/rules/SunOS64Gcc/mplib        |  3 +++
 wmake/rules/SunOS64Gcc/mplibFJMPI   |  3 +++
 wmake/rules/SunOS64Gcc/mplibOPENMPI |  3 +++
 18 files changed, 119 insertions(+), 2 deletions(-)
 create mode 100644 wmake/rules/SunOS64Gcc/X
 create mode 100644 wmake/rules/SunOS64Gcc/c
 create mode 100644 wmake/rules/SunOS64Gcc/c++
 create mode 100644 wmake/rules/SunOS64Gcc/c++Debug
 create mode 100644 wmake/rules/SunOS64Gcc/c++Opt
 create mode 100644 wmake/rules/SunOS64Gcc/c++Prof
 create mode 100644 wmake/rules/SunOS64Gcc/cDebug
 create mode 100644 wmake/rules/SunOS64Gcc/cOpt
 create mode 100644 wmake/rules/SunOS64Gcc/cProf
 create mode 100644 wmake/rules/SunOS64Gcc/general
 create mode 100644 wmake/rules/SunOS64Gcc/mplib
 create mode 100644 wmake/rules/SunOS64Gcc/mplibFJMPI
 create mode 100644 wmake/rules/SunOS64Gcc/mplibOPENMPI

diff --git a/etc/bashrc b/etc/bashrc
index 7c9f9ecad3..e9f9c6f9cf 100644
--- a/etc/bashrc
+++ b/etc/bashrc
@@ -172,6 +172,17 @@ Linux)
     esac
     ;;
 
+SunOS)
+    WM_ARCH=SunOS64
+    export WM_COMPILER_LIB_ARCH=64
+    export WM_CC='gcc'
+    export WM_CXX='g++'
+    export WM_CFLAGS='-mabi=64 -fPIC'
+    export WM_CXXFLAGS='-mabi=64 -fPIC'
+    export WM_LDFLAGS='-mabi=64 -G0'
+    export WM_MPLIB=FJMPI
+    ;;
+
 *)    	# an unsupported operating system
     cat <<USAGE
 
diff --git a/etc/cshrc b/etc/cshrc
index b1dfd9ec34..5952f6d3c5 100644
--- a/etc/cshrc
+++ b/etc/cshrc
@@ -149,7 +149,7 @@ case Linux:
         setenv WM_ARCH linuxIA64
         setenv WM_COMPILER I64
         breaksw
-    mips64)
+    case mips64:
         setenv WM_ARCH SiCortex64
         setenv WM_COMPILER_LIB_ARCH 64
         setenv WM_CC 'gcc'
@@ -165,6 +165,17 @@ case Linux:
     endsw
     breaksw
 
+case SunOS:
+    setenv WM_ARCH SunOS64
+    setenv WM_COMPILER_LIB_ARCH 64
+    setenv WM_CC 'gcc'
+    setenv WM_CXX 'g++'
+    setenv WM_CFLAGS '-mabi=64 -fPIC'
+    setenv WM_CXXFLAGS '-mabi=64 -fPIC'
+    setenv WM_LDFLAGS '-mabi=64 -G0'
+    setenv WM_MPLIB FJMPI
+    breaksw
+
 default:
     echo
     echo "Your '$WM_ARCH' operating system is not supported by this release"
diff --git a/etc/settings.csh b/etc/settings.csh
index 73b14eb38c..a38ec5eacd 100644
--- a/etc/settings.csh
+++ b/etc/settings.csh
@@ -222,6 +222,15 @@ case MPI:
     setenv FOAM_MPI_LIBBIN $FOAM_LIBBIN/mpi
     breaksw
 
+case FJMPI:
+    setenv MPI_ARCH_PATH /opt/FJSVmpi2
+    setenv FOAM_MPI_LIBBIN $FOAM_LIBBIN/mpi
+    _foamAddPath $MPI_ARCH_PATH/bin
+    _foamAddLib  $MPI_ARCH_PATH/lib/sparcv9
+    _foamAddLib  /opt/FSUNf90/lib/sparcv9
+    _foamAddLib  /opt/FJSVpnidt/lib
+    breaksw
+
 default:
     setenv FOAM_MPI_LIBBIN $FOAM_LIBBIN/dummy
     breaksw
diff --git a/etc/settings.sh b/etc/settings.sh
index ff56af0b95..7df3f6d590 100644
--- a/etc/settings.sh
+++ b/etc/settings.sh
@@ -249,6 +249,16 @@ MPI)
     export FOAM_MPI_LIBBIN=$FOAM_LIBBIN/mpi
     ;;
 
+FJMPI)
+    export MPI_ARCH_PATH=/opt/FJSVmpi2
+    export FOAM_MPI_LIBBIN=$FOAM_LIBBIN/mpi
+
+    _foamAddPath $MPI_ARCH_PATH/bin
+    _foamAddLib  $MPI_ARCH_PATH/lib/sparcv9
+    _foamAddLib  /opt/FSUNf90/lib/sparcv9
+    _foamAddLib  /opt/FJSVpnidt/lib
+    ;;
+
 *)
     export FOAM_MPI_LIBBIN=$FOAM_LIBBIN/dummy
     ;;
diff --git a/src/OSspecific/Unix/Make/files b/src/OSspecific/Unix/Make/files
index f83513ac4a..b6e32d80d6 100644
--- a/src/OSspecific/Unix/Make/files
+++ b/src/OSspecific/Unix/Make/files
@@ -8,7 +8,11 @@ fileStat.C
 Unix.C
 cpuTime/cpuTime.C
 clockTime/clockTime.C
+
+#ifndef SunOS64
 printStack.C
-/*dummyPrintStack.C*/
+#else
+dummyPrintStack.C
+#endif
 
 LIB = $(FOAM_LIBBIN)/libOSspecific
diff --git a/wmake/rules/SunOS64Gcc/X b/wmake/rules/SunOS64Gcc/X
new file mode 100644
index 0000000000..5d1f9c5cc5
--- /dev/null
+++ b/wmake/rules/SunOS64Gcc/X
@@ -0,0 +1,3 @@
+XFLAGS     =
+XINC       = $(XFLAGS) -I/usr/X11R6/include
+XLIBS      = -L/usr/X11R6/lib64 -lXext -lX11
diff --git a/wmake/rules/SunOS64Gcc/c b/wmake/rules/SunOS64Gcc/c
new file mode 100644
index 0000000000..80bb80f32f
--- /dev/null
+++ b/wmake/rules/SunOS64Gcc/c
@@ -0,0 +1,16 @@
+.SUFFIXES: .c .h
+
+cWARN        = -Wall
+
+cc          = gcc -m64
+
+include $(RULES)/c$(WM_COMPILE_OPTION)
+
+cFLAGS      = $(GFLAGS) $(cWARN) $(cOPT) $(cDBUG) $(LIB_HEADER_DIRS) -fPIC
+
+ctoo        = $(WM_SCHEDULER) $(cc) $(cFLAGS) -c $$SOURCE -o $@
+
+LINK_LIBS   = $(cDBUG)
+
+LINKLIBSO   = $(cc) -shared
+LINKEXE     = $(cc) -Xlinker -z -Xlinker nodefs
diff --git a/wmake/rules/SunOS64Gcc/c++ b/wmake/rules/SunOS64Gcc/c++
new file mode 100644
index 0000000000..36b5c29b20
--- /dev/null
+++ b/wmake/rules/SunOS64Gcc/c++
@@ -0,0 +1,21 @@
+.SUFFIXES: .C .cxx .cc .cpp
+
+c++WARN     = -Wall -Wno-strict-aliasing -Wextra -Wno-unused-parameter -Wold-style-cast
+
+CC          = g++ -m64
+
+include $(RULES)/c++$(WM_COMPILE_OPTION)
+
+ptFLAGS     = -DNoRepository -ftemplate-depth-40
+
+c++FLAGS    = $(GFLAGS) $(c++WARN) $(c++OPT) $(c++DBUG) $(ptFLAGS) $(LIB_HEADER_DIRS) -fPIC
+
+Ctoo        = $(WM_SCHEDULER) $(CC) $(c++FLAGS) -c $$SOURCE -o $@
+cxxtoo      = $(Ctoo)
+cctoo       = $(Ctoo)
+cpptoo      = $(Ctoo)
+
+LINK_LIBS   = $(c++DBUG)
+
+LINKLIBSO   = $(CC) $(c++FLAGS) -shared
+LINKEXE     = $(CC) $(c++FLAGS)
diff --git a/wmake/rules/SunOS64Gcc/c++Debug b/wmake/rules/SunOS64Gcc/c++Debug
new file mode 100644
index 0000000000..19bdb9c334
--- /dev/null
+++ b/wmake/rules/SunOS64Gcc/c++Debug
@@ -0,0 +1,2 @@
+c++DBUG    = -ggdb3 -DFULLDEBUG
+c++OPT      = -O0 -fdefault-inline
diff --git a/wmake/rules/SunOS64Gcc/c++Opt b/wmake/rules/SunOS64Gcc/c++Opt
new file mode 100644
index 0000000000..4cb95b5b3a
--- /dev/null
+++ b/wmake/rules/SunOS64Gcc/c++Opt
@@ -0,0 +1,2 @@
+c++DBUG     = 
+c++OPT      = -O3
diff --git a/wmake/rules/SunOS64Gcc/c++Prof b/wmake/rules/SunOS64Gcc/c++Prof
new file mode 100644
index 0000000000..3bda4dad55
--- /dev/null
+++ b/wmake/rules/SunOS64Gcc/c++Prof
@@ -0,0 +1,2 @@
+c++DBUG    = -pg
+c++OPT     = -O2
diff --git a/wmake/rules/SunOS64Gcc/cDebug b/wmake/rules/SunOS64Gcc/cDebug
new file mode 100644
index 0000000000..72b638f458
--- /dev/null
+++ b/wmake/rules/SunOS64Gcc/cDebug
@@ -0,0 +1,2 @@
+cDBUG       = -ggdb -DFULLDEBUG
+cOPT        = -O1 -fdefault-inline -finline-functions
diff --git a/wmake/rules/SunOS64Gcc/cOpt b/wmake/rules/SunOS64Gcc/cOpt
new file mode 100644
index 0000000000..aaaebef3d3
--- /dev/null
+++ b/wmake/rules/SunOS64Gcc/cOpt
@@ -0,0 +1,2 @@
+cDBUG       = 
+cOPT        = -O3 -fno-gcse
diff --git a/wmake/rules/SunOS64Gcc/cProf b/wmake/rules/SunOS64Gcc/cProf
new file mode 100644
index 0000000000..ca3ac9bf5f
--- /dev/null
+++ b/wmake/rules/SunOS64Gcc/cProf
@@ -0,0 +1,2 @@
+cDBUG       = -pg
+cOPT        = -O2
diff --git a/wmake/rules/SunOS64Gcc/general b/wmake/rules/SunOS64Gcc/general
new file mode 100644
index 0000000000..70ee588d82
--- /dev/null
+++ b/wmake/rules/SunOS64Gcc/general
@@ -0,0 +1,11 @@
+CPP        = /lib/cpp $(GFLAGS)
+LD         = ld -64
+
+PROJECT_LIBS = -l$(WM_PROJECT) -liberty -lnsl -lsocket -L$(FOAM_LIBBIN)/dummy -lPstream 
+
+include $(GENERAL_RULES)/standard
+
+include $(RULES)/X
+include $(RULES)/c
+include $(RULES)/c++
+include $(GENERAL_RULES)/cint
diff --git a/wmake/rules/SunOS64Gcc/mplib b/wmake/rules/SunOS64Gcc/mplib
new file mode 100644
index 0000000000..8a84b40146
--- /dev/null
+++ b/wmake/rules/SunOS64Gcc/mplib
@@ -0,0 +1,3 @@
+PFLAGS     = 
+PINC       = 
+PLIBS      = 
diff --git a/wmake/rules/SunOS64Gcc/mplibFJMPI b/wmake/rules/SunOS64Gcc/mplibFJMPI
new file mode 100644
index 0000000000..2d6c16797a
--- /dev/null
+++ b/wmake/rules/SunOS64Gcc/mplibFJMPI
@@ -0,0 +1,3 @@
+PFLAGS     =
+PINC       = -I$(MPI_ARCH_PATH)/include/sparcv9/mpi -DMPIPP_H
+PLIBS      = -L$(MPI_ARCH_PATH)/lib/sparcv9 -lmpi_f -L/opt/FSUNf90/lib/sparcv9 -lfjgmp64 -L/opt/FJSVpnidt/lib -lfjidt -ljrm -lfj90i -lfj90f -lfj90fmt  -lelf
diff --git a/wmake/rules/SunOS64Gcc/mplibOPENMPI b/wmake/rules/SunOS64Gcc/mplibOPENMPI
new file mode 100644
index 0000000000..834d2d3e22
--- /dev/null
+++ b/wmake/rules/SunOS64Gcc/mplibOPENMPI
@@ -0,0 +1,3 @@
+PFLAGS     = -DOMPI_SKIP_MPICXX
+PINC       = -I$(MPI_ARCH_PATH)/include
+PLIBS      = -L$(MPI_ARCH_PATH)/lib -lmpi

From 8fa4620e0afe905fe1d896ac5ee70d54476e8abd Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Mon, 2 Feb 2009 12:24:32 +0000
Subject: [PATCH 5/5] cutting plane

---
 .../postProcessing/sampling/sample/sampleDict |  19 +
 src/sampling/Make/files                       |   1 +
 .../sampledCuttingPlane/sampledCuttingPlane.C | 420 ++++++++++++++++++
 .../sampledCuttingPlane/sampledCuttingPlane.H | 258 +++++++++++
 .../sampledCuttingPlaneTemplates.C            |  97 ++++
 5 files changed, 795 insertions(+)
 create mode 100644 src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
 create mode 100644 src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
 create mode 100644 src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C

diff --git a/applications/utilities/postProcessing/sampling/sample/sampleDict b/applications/utilities/postProcessing/sampling/sample/sampleDict
index 1616bbb662..b605f29e57 100644
--- a/applications/utilities/postProcessing/sampling/sample/sampleDict
+++ b/applications/utilities/postProcessing/sampling/sample/sampleDict
@@ -177,6 +177,25 @@ surfaces
         interpolate     false;
         regularise      false;              // do not simplify
     }
+
+    triangleCut
+    {
+        // Cutingplaneusing iso surface
+        type            cuttingPlane;
+        planeType       pointAndNormal;
+        pointAndNormalDict
+        {
+            basePoint       (0.4 0 0.4);
+            normalVector    (1 0.2 0.2);
+        }
+        interpolate     true;
+
+        //zone            ABC;          // Optional: zone only
+        //exposedPatchName fixedWalls;  // Optional: zone only
+
+        // regularise      false;    // Optional: do not simplify
+    }
+
 );
 
 
diff --git a/src/sampling/Make/files b/src/sampling/Make/files
index 43bc910a3c..004f81d435 100644
--- a/src/sampling/Make/files
+++ b/src/sampling/Make/files
@@ -29,6 +29,7 @@ sampledSurface/isoSurface/sampledIsoSurface.C
 sampledSurface/isoSurface/isoSurfaceCell.C
 sampledSurface/isoSurface/sampledIsoSurfaceCell.C
 sampledSurface/distanceSurface/distanceSurface.C
+sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
 sampledSurface/sampledSurface/sampledSurface.C
 sampledSurface/sampledSurfaces/sampledSurfaces.C
 sampledSurface/sampledSurfacesFunctionObject/sampledSurfacesFunctionObject.C
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
new file mode 100644
index 0000000000..85dbeca595
--- /dev/null
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
@@ -0,0 +1,420 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "sampledCuttingPlane.H"
+#include "dictionary.H"
+#include "volFields.H"
+#include "volPointInterpolation.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvMesh.H"
+#include "isoSurface.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(sampledCuttingPlane, 0);
+    addNamedToRunTimeSelectionTable(sampledSurface, sampledCuttingPlane, word, cuttingPlane);
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::sampledCuttingPlane::createGeometry()
+{
+    if (debug)
+    {
+        Pout<< "sampledCuttingPlane::createGeometry :updating geometry."
+            << endl;
+    }
+
+    // Clear any stored topologies
+    facesPtr_.clear();
+    isoSurfPtr_.ptr();
+    pointDistance_.clear();
+    cellDistancePtr_.clear();
+
+
+    // Get any subMesh
+    if (zoneID_.index() != -1 && !subMeshPtr_.valid())
+    {
+        const polyBoundaryMesh& patches = mesh().boundaryMesh();
+
+        // Patch to put exposed internal faces into
+        label exposedPatchI = patches.findPatchID(exposedPatchName_);
+
+        if (debug)
+        {
+            Info<< "Allocating subset of size "
+                << mesh().cellZones()[zoneID_.index()].size()
+                << " with exposed faces into patch "
+                << patches[exposedPatchI].name() << endl;
+        }
+
+        subMeshPtr_.reset
+        (
+            new fvMeshSubset(static_cast<const fvMesh&>(mesh()))
+        );
+        subMeshPtr_().setLargeCellSubset
+        (
+            labelHashSet(mesh().cellZones()[zoneID_.index()]),
+            exposedPatchI
+        );
+    }
+
+
+    // Select either the submesh or the underlying mesh
+    const fvMesh& fvm =
+    (
+        subMeshPtr_.valid()
+      ? subMeshPtr_().subMesh()
+      : static_cast<const fvMesh&>(mesh())
+    );
+
+
+    // Distance to cell centres
+    // ~~~~~~~~~~~~~~~~~~~~~~~~
+
+    cellDistancePtr_.reset
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "cellDistance",
+                fvm.time().timeName(),
+                fvm.time(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            fvm,
+            dimensionedScalar("zero", dimLength, 0)
+        )
+    );
+    volScalarField& cellDistance = cellDistancePtr_();
+
+    // Internal field
+    {
+        const pointField& cc = fvm.C();
+        scalarField& fld = cellDistance.internalField();
+
+        forAll(cc, i)
+        {
+            // Signed distance
+            fld[i] = (cc[i] - plane_.refPoint()) & plane_.normal();
+        }
+    }
+
+    // Patch fields
+    {
+        forAll(fvm.C().boundaryField(), patchI)
+        {
+            const pointField& cc = fvm.C().boundaryField()[patchI];
+            fvPatchScalarField& fld = cellDistance.boundaryField()[patchI];
+
+            forAll(fld, i)
+            {
+                fld[i] =  (cc[i] - plane_.refPoint()) & plane_.normal();
+            }
+        }
+    }
+
+
+    // On processor patches the mesh.C() will already be the cell centre
+    // on the opposite side so no need to swap cellDistance.
+
+
+    // Distance to points
+    pointDistance_.setSize(fvm.nPoints());
+    {
+        const pointField& pts = fvm.points();
+
+        forAll(pointDistance_, i)
+        {
+            pointDistance_[i] = (pts[i] - plane_.refPoint()) & plane_.normal();
+        }
+    }
+
+
+    if (debug)
+    {
+        Pout<< "Writing cell distance:" << cellDistance.objectPath() << endl;
+        cellDistance.write();
+        pointScalarField pDist
+        (
+            IOobject
+            (
+                "pointDistance",
+                fvm.time().timeName(),
+                fvm.time(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            pointMesh::New(fvm),
+            dimensionedScalar("zero", dimLength, 0)
+        );
+        pDist.internalField() = pointDistance_;
+
+        Pout<< "Writing point distance:" << pDist.objectPath() << endl;
+        pDist.write();
+    }
+
+
+    //- Direct from cell field and point field.
+    isoSurfPtr_.reset
+    (
+        new isoSurface
+        (
+            cellDistance,
+            pointDistance_,
+            0.0,
+            regularise_
+        )
+    );
+
+    if (debug)
+    {
+        print(Pout);
+        Pout<< endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::sampledCuttingPlane::sampledCuttingPlane
+(
+    const word& name,
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    sampledSurface(name, mesh, dict),
+    plane_(dict),
+    mergeTol_(dict.lookupOrDefault("mergeTol", 1E-6)),
+    regularise_(dict.lookupOrDefault("regularise", true)),
+    zoneID_(dict.lookupOrDefault("zone", word::null), mesh.cellZones()),
+    exposedPatchName_(word::null),
+    needsUpdate_(true),
+    subMeshPtr_(NULL),
+    cellDistancePtr_(NULL),
+    isoSurfPtr_(NULL),
+    facesPtr_(NULL)
+{
+    if (zoneID_.index() != -1)
+    {
+        dict.lookup("exposedPatchName") >> exposedPatchName_;
+
+        if (mesh.boundaryMesh().findPatchID(exposedPatchName_) == -1)
+        {
+            FatalErrorIn
+            (
+                "sampledCuttingPlane::sampledCuttingPlane"
+                "(const word&, const polyMesh&, const dictionary&)"
+            )   << "Cannot find patch " << exposedPatchName_
+                << " in which to put exposed faces." << endl
+                << "Valid patches are " << mesh.boundaryMesh().names()
+                << exit(FatalError);
+        }
+
+        if (debug && zoneID_.index() != -1)
+        {
+            Info<< "Restricting to cellZone " << zoneID_.name()
+                << " with exposed internal faces into patch "
+                << exposedPatchName_ << endl;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::sampledCuttingPlane::~sampledCuttingPlane()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::sampledCuttingPlane::needsUpdate() const
+{
+    return needsUpdate_;
+}
+
+
+bool Foam::sampledCuttingPlane::expire()
+{
+    if (debug)
+    {
+        Pout<< "sampledCuttingPlane::expire :"
+            << " have-facesPtr_:" << facesPtr_.valid()
+            << " needsUpdate_:" << needsUpdate_ << endl;
+    }
+
+    // Clear any stored topologies
+    facesPtr_.clear();
+
+    // already marked as expired
+    if (needsUpdate_)
+    {
+        return false;
+    }
+
+    needsUpdate_ = true;
+    return true;
+}
+
+
+bool Foam::sampledCuttingPlane::update()
+{
+    if (debug)
+    {
+        Pout<< "sampledCuttingPlane::update :"
+            << " have-facesPtr_:" << facesPtr_.valid()
+            << " needsUpdate_:" << needsUpdate_ << endl;
+    }
+
+    if (!needsUpdate_)
+    {
+        return false;
+    }
+
+    createGeometry();
+
+    needsUpdate_ = false;
+    return true;
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::sampledCuttingPlane::sample
+(
+    const volScalarField& vField
+) const
+{
+    return sampleField(vField);
+}
+
+
+Foam::tmp<Foam::vectorField>
+Foam::sampledCuttingPlane::sample
+(
+    const volVectorField& vField
+) const
+{
+    return sampleField(vField);
+}
+
+
+Foam::tmp<Foam::sphericalTensorField>
+Foam::sampledCuttingPlane::sample
+(
+    const volSphericalTensorField& vField
+) const
+{
+    return sampleField(vField);
+}
+
+
+Foam::tmp<Foam::symmTensorField>
+Foam::sampledCuttingPlane::sample
+(
+    const volSymmTensorField& vField
+) const
+{
+    return sampleField(vField);
+}
+
+
+Foam::tmp<Foam::tensorField>
+Foam::sampledCuttingPlane::sample
+(
+    const volTensorField& vField
+) const
+{
+    return sampleField(vField);
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::sampledCuttingPlane::interpolate
+(
+    const interpolation<scalar>& interpolator
+) const
+{
+    return interpolateField(interpolator);
+}
+
+
+Foam::tmp<Foam::vectorField>
+Foam::sampledCuttingPlane::interpolate
+(
+    const interpolation<vector>& interpolator
+) const
+{
+    return interpolateField(interpolator);
+}
+
+Foam::tmp<Foam::sphericalTensorField>
+Foam::sampledCuttingPlane::interpolate
+(
+    const interpolation<sphericalTensor>& interpolator
+) const
+{
+    return interpolateField(interpolator);
+}
+
+
+Foam::tmp<Foam::symmTensorField>
+Foam::sampledCuttingPlane::interpolate
+(
+    const interpolation<symmTensor>& interpolator
+) const
+{
+    return interpolateField(interpolator);
+}
+
+
+Foam::tmp<Foam::tensorField>
+Foam::sampledCuttingPlane::interpolate
+(
+    const interpolation<tensor>& interpolator
+) const
+{
+    return interpolateField(interpolator);
+}
+
+
+void Foam::sampledCuttingPlane::print(Ostream& os) const
+{
+    os  << "sampledCuttingPlane: " << name() << " :"
+        << "  plane:" << plane_
+        << "  faces:" << faces().size()
+        << "  points:" << points().size();
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
new file mode 100644
index 0000000000..fcf14db8bd
--- /dev/null
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
@@ -0,0 +1,258 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::sampledCuttingPlane
+
+Description
+    A sampledSurface defined by a plane
+
+SourceFiles
+    sampledCuttingPlane.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef sampledCuttingPlane_H
+#define sampledCuttingPlane_H
+
+#include "sampledSurface.H"
+#include "isoSurface.H"
+#include "plane.H"
+#include "ZoneIDs.H"
+#include "fvMeshSubset.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class sampledCuttingPlane Declaration
+\*---------------------------------------------------------------------------*/
+
+class sampledCuttingPlane
+:
+    public sampledSurface
+{
+    // Private data
+
+        //- Plane
+        const plane plane_;
+
+        //- Merge tolerance
+        const scalar mergeTol_;
+
+        //- Whether to coarsen
+        const Switch regularise_;
+
+        //- zone name/index (if restricted to zones)
+        mutable cellZoneID zoneID_;
+
+        //- for zones: patch to put exposed faces into
+        mutable word exposedPatchName_;
+
+        //- Track if the surface needs an update
+        mutable bool needsUpdate_;
+
+
+        //- Optional subsetted mesh
+        autoPtr<fvMeshSubset> subMeshPtr_;
+
+        //- Distance to cell centres
+        autoPtr<volScalarField> cellDistancePtr_;
+
+        //- Distance to points
+        scalarField pointDistance_;
+
+        //- Constructed iso surface
+        autoPtr<isoSurface> isoSurfPtr_;
+
+        //- triangles converted to faceList
+        mutable autoPtr<faceList> facesPtr_;
+
+
+    // Private Member Functions
+
+        //- Create iso surface
+        void createGeometry();
+
+        //- sample field on faces
+        template <class Type>
+        tmp<Field<Type> > sampleField
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& vField
+        ) const;
+
+
+        template <class Type>
+        tmp<Field<Type> >
+        interpolateField(const interpolation<Type>&) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("sampledCuttingPlane");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        sampledCuttingPlane
+        (
+            const word& name,
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+
+    // Destructor
+
+        virtual ~sampledCuttingPlane();
+
+
+    // Member Functions
+
+        //- Does the surface need an update?
+        virtual bool needsUpdate() const;
+
+        //- Mark the surface as needing an update.
+        //  May also free up unneeded data.
+        //  Return false if surface was already marked as expired.
+        virtual bool expire();
+
+        //- Update the surface as required.
+        //  Do nothing (and return false) if no update was needed
+        virtual bool update();
+
+        //- Points of surface
+        virtual const pointField& points() const
+        {
+            return surface().points();
+        }
+
+        //- Faces of surface
+        virtual const faceList& faces() const
+        {
+            if (facesPtr_.empty())
+            {
+                const triSurface& s = surface();
+
+                facesPtr_.reset(new faceList(s.size()));
+
+                forAll(s, i)
+                {
+                    facesPtr_()[i] = s[i].triFaceFace();
+                }
+            }
+            return facesPtr_;
+        }
+
+
+        const isoSurface& surface() const
+        {
+            return isoSurfPtr_();
+        }
+
+        //- sample field on surface
+        virtual tmp<scalarField> sample
+        (
+            const volScalarField&
+        ) const;
+
+        //- sample field on surface
+        virtual tmp<vectorField> sample
+        (
+            const volVectorField&
+        ) const;
+
+        //- sample field on surface
+        virtual tmp<sphericalTensorField> sample
+        (
+            const volSphericalTensorField&
+        ) const;
+
+        //- sample field on surface
+        virtual tmp<symmTensorField> sample
+        (
+            const volSymmTensorField&
+        ) const;
+
+        //- sample field on surface
+        virtual tmp<tensorField> sample
+        (
+            const volTensorField&
+        ) const;
+
+
+        //- interpolate field on surface
+        virtual tmp<scalarField> interpolate
+        (
+            const interpolation<scalar>&
+        ) const;
+
+        //- interpolate field on surface
+        virtual tmp<vectorField> interpolate
+        (
+            const interpolation<vector>&
+        ) const;
+
+        //- interpolate field on surface
+        virtual tmp<sphericalTensorField> interpolate
+        (
+            const interpolation<sphericalTensor>&
+        ) const;
+
+        //- interpolate field on surface
+        virtual tmp<symmTensorField> interpolate
+        (
+            const interpolation<symmTensor>&
+        ) const;
+
+        //- interpolate field on surface
+        virtual tmp<tensorField> interpolate
+        (
+            const interpolation<tensor>&
+        ) const;
+
+        //- Write
+        virtual void print(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "sampledCuttingPlaneTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C
new file mode 100644
index 0000000000..f83aeca4e5
--- /dev/null
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C
@@ -0,0 +1,97 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "volPointInterpolation.H"
+#include "sampledCuttingPlane.H"
+#include "isoSurface.H"
+#include "volFieldsFwd.H"
+#include "pointFields.H"
+#include "volPointInterpolation.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template <class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::sampledCuttingPlane::sampleField
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vField
+) const
+{
+    return tmp<Field<Type> >(new Field<Type>(vField, surface().meshCells()));
+}
+
+
+template <class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::sampledCuttingPlane::interpolateField
+(
+    const interpolation<Type>& interpolator
+) const
+{
+    // Get fields to sample. Assume volPointInterpolation!
+    const GeometricField<Type, fvPatchField, volMesh>& volFld =
+        interpolator.psi();
+
+    if (subMeshPtr_.valid())
+    {
+        tmp<GeometricField<Type, fvPatchField, volMesh> > tvolSubFld =
+            subMeshPtr_().interpolate(volFld);
+
+        const GeometricField<Type, fvPatchField, volMesh>& volSubFld =
+            tvolSubFld();
+
+        tmp<GeometricField<Type, pointPatchField, pointMesh> > tpointSubFld =
+            volPointInterpolation::New(volSubFld.mesh()).interpolate(volSubFld);
+
+        // Sample.
+        return surface().interpolate
+        (
+            cellDistancePtr_(),
+            pointDistance_,
+            volSubFld,
+            tpointSubFld()
+        );
+    }
+    else
+    {
+        tmp<GeometricField<Type, pointPatchField, pointMesh> > tpointFld
+        (
+            volPointInterpolation::New(volFld.mesh()).interpolate(volFld)
+        );
+
+        // Sample.
+        return surface().interpolate
+        (
+            cellDistancePtr_(),
+            pointDistance_,
+            volFld,
+            tpointFld()
+        );
+    }
+}
+
+
+// ************************************************************************* //