diff --git a/applications/test/BSpline/BSplineTest.C b/applications/test/BSpline/BSplineTest.C
index cb858d5560..91f3efad23 100644
--- a/applications/test/BSpline/BSplineTest.C
+++ b/applications/test/BSpline/BSplineTest.C
@@ -25,7 +25,6 @@ License
 \*---------------------------------------------------------------------------*/
 #include "argList.H"
 
-#include "simpleMatrix.H"
 #include "vector.H"
 #include "IFstream.H"
 #include "BSpline.H"
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkMesh.C b/applications/utilities/mesh/manipulation/checkMesh/checkMesh.C
index 96fb76ffe2..6caeffbd4e 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkMesh.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkMesh.C
@@ -94,7 +94,8 @@ int main(int argc, char *argv[])
 
             noFailedChecks += checkGeometry(mesh, allGeometry);
 
-            reduce(noFailedChecks, sumOp<label>());
+            // Note: no reduction in noFailedChecks necessary since is
+            //       counter of checks, not counter of failed cells,faces etc.
 
             if (noFailedChecks == 0)
             {
@@ -112,8 +113,6 @@ int main(int argc, char *argv[])
 
             label nFailedChecks = checkGeometry(mesh, allGeometry);
 
-            reduce(nFailedChecks, sumOp<label>());
-
             if (nFailedChecks)
             {
                 Info<< "\nFailed " << nFailedChecks << " mesh checks.\n"
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/writeFuns.C b/applications/utilities/postProcessing/dataConversion/foamToVTK/writeFuns.C
index f699e17588..82a8df0cca 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/writeFuns.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/writeFuns.C
@@ -220,7 +220,7 @@ void Foam::writeFuns::writePointDataHeader
 }
 
 
-void Foam::writeFuns::insert(const scalar& pt, DynamicList<floatScalar>& dest)
+void Foam::writeFuns::insert(const scalar pt, DynamicList<floatScalar>& dest)
 {
     dest.append(float(pt));
 }
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/writeFuns.H b/applications/utilities/postProcessing/dataConversion/foamToVTK/writeFuns.H
index cf7a4e992c..c38bbda669 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/writeFuns.H
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/writeFuns.H
@@ -86,7 +86,7 @@ public:
 
     // Convert to VTK and store
 
-        static void insert(const scalar&, DynamicList<floatScalar>&);
+        static void insert(const scalar, DynamicList<floatScalar>&);
         static void insert(const point&, DynamicList<floatScalar>&);
         static void insert(const sphericalTensor&, DynamicList<floatScalar>&);
         static void insert(const symmTensor&, DynamicList<floatScalar>&);
diff --git a/src/OpenFOAM/db/scalarRange/scalarRange.C b/src/OpenFOAM/db/scalarRange/scalarRange.C
index 7623e14c27..393b19bba4 100644
--- a/src/OpenFOAM/db/scalarRange/scalarRange.C
+++ b/src/OpenFOAM/db/scalarRange/scalarRange.C
@@ -42,7 +42,7 @@ Foam::scalarRange::scalarRange()
 {}
 
 
-Foam::scalarRange::scalarRange(const scalar& lower, const scalar& upper)
+Foam::scalarRange::scalarRange(const scalar lower, const scalar upper)
 :
     type_(RANGE),
     value_(lower),
@@ -123,7 +123,7 @@ Foam::scalar Foam::scalarRange::upper() const
 }
 
 
-bool Foam::scalarRange::selected(const scalar& value) const
+bool Foam::scalarRange::selected(const scalar value) const
 {
     switch (type_)
     {
diff --git a/src/OpenFOAM/db/scalarRange/scalarRange.H b/src/OpenFOAM/db/scalarRange/scalarRange.H
index b996871664..f123809237 100644
--- a/src/OpenFOAM/db/scalarRange/scalarRange.H
+++ b/src/OpenFOAM/db/scalarRange/scalarRange.H
@@ -92,7 +92,7 @@ public:
         scalarRange();
 
         //- Construct a Range
-        scalarRange(const scalar& lower, const scalar& upper);
+        scalarRange(const scalar lower, const scalar upper);
 
         //- Construct from Istream.
         //  Since commas can be used as list delimiters,
@@ -119,7 +119,7 @@ public:
         scalar upper() const;
 
         //- Return true if the value is within the range
-        bool selected(const scalar&) const;
+        bool selected(const scalar) const;
 
 
     // Member Operators
diff --git a/src/OpenFOAM/db/scalarRange/scalarRanges.C b/src/OpenFOAM/db/scalarRange/scalarRanges.C
index c59e2b5a9d..f8b3c077b1 100644
--- a/src/OpenFOAM/db/scalarRange/scalarRanges.C
+++ b/src/OpenFOAM/db/scalarRange/scalarRanges.C
@@ -57,7 +57,7 @@ Foam::scalarRanges::scalarRanges(Istream& is)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-bool Foam::scalarRanges::selected(const scalar& value) const
+bool Foam::scalarRanges::selected(const scalar value) const
 {
     forAll(*this, i)
     {
diff --git a/src/OpenFOAM/db/scalarRange/scalarRanges.H b/src/OpenFOAM/db/scalarRange/scalarRanges.H
index 5483454ae6..c3bf1ea668 100644
--- a/src/OpenFOAM/db/scalarRange/scalarRanges.H
+++ b/src/OpenFOAM/db/scalarRange/scalarRanges.H
@@ -66,7 +66,7 @@ public:
     // Member Functions
 
         //- Return true if the given value is within the ranges
-        bool selected(const scalar&) const;
+        bool selected(const scalar) const;
 
         //- Return the set of selected entries in the given list
         //  that are within the ranges
diff --git a/src/OpenFOAM/global/unitConversion/unitConversion.H b/src/OpenFOAM/global/unitConversion/unitConversion.H
index 73d80dd57e..c83582cf84 100644
--- a/src/OpenFOAM/global/unitConversion/unitConversion.H
+++ b/src/OpenFOAM/global/unitConversion/unitConversion.H
@@ -43,13 +43,13 @@ namespace Foam
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 //- Conversion from degrees to radians
-inline scalar degToRad(const scalar& deg)
+inline scalar degToRad(const scalar deg)
 {
     return (deg * constant::mathematical::pi/180.0);
 }
 
 //- Conversion from radians to degrees
-inline scalar radToDeg(const scalar& rad)
+inline scalar radToDeg(const scalar rad)
 {
     return (rad * 180.0/constant::mathematical::pi);
 }
diff --git a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C
index 416f06c1bf..f412f20ee0 100644
--- a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C
+++ b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C
@@ -31,16 +31,21 @@ License
 template<class Type>
 Foam::simpleMatrix<Type>::simpleMatrix(const label mSize)
 :
-    scalarSquareMatrix(mSize, mSize, pTraits<scalar>::zero),
-    source_(mSize, pTraits<Type>::zero)
+    scalarSquareMatrix(mSize),
+    source_(mSize)
 {}
 
 
 template<class Type>
-Foam::simpleMatrix<Type>::simpleMatrix(const label mSize, const scalar& t)
+Foam::simpleMatrix<Type>::simpleMatrix
+(
+    const label mSize,
+    const scalar coeffVal,
+    const Type& sourceVal
+)
 :
-    scalarSquareMatrix(mSize, mSize, t),
-    source_(mSize, pTraits<Type>::zero)
+    scalarSquareMatrix(mSize, mSize, coeffVal),
+    source_(mSize, sourceVal)
 {}
 
 
diff --git a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H
index e203972ffe..a32f318d33 100644
--- a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H
+++ b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H
@@ -74,11 +74,12 @@ public:
 
     // Constructors
 
-        //- Construct given size, initializing matrix coefficients to zero
+        //- Construct given size
+        //  Note: this does not initialise the coefficients or the source.
         simpleMatrix(const label);
 
-        //- Construct given size, and initial value for matrix coefficients
-        simpleMatrix(const label, const scalar&);
+        //- Construct given size and initial values for coefficients and source
+        simpleMatrix(const label, const scalar, const Type&);
 
         //- Construct from components
         simpleMatrix(const scalarSquareMatrix&, const Field<Type>&);
diff --git a/src/OpenFOAM/matrices/solution/solution.C b/src/OpenFOAM/matrices/solution/solution.C
index c63dceec4e..ba0db67bd4 100644
--- a/src/OpenFOAM/matrices/solution/solution.C
+++ b/src/OpenFOAM/matrices/solution/solution.C
@@ -60,6 +60,7 @@ Foam::solution::solution(const objectRegistry& obr, const fileName& dictName)
         )
     ),
     cache_(ITstream("cache", tokenList())()),
+    caching_(false),
     relaxationFactors_(ITstream("relaxationFactors", tokenList())()),
     defaultRelaxationFactor_(0),
     solvers_(ITstream("solvers", tokenList())())
@@ -150,12 +151,19 @@ Foam::label Foam::solution::upgradeSolverDict
 
 bool Foam::solution::cache(const word& name) const
 {
-    if (debug)
+    if (caching_)
     {
-        Info<< "Find cache entry for " << name << endl;
-    }
+        if (debug)
+        {
+            Info<< "Cache: find entry for " << name << endl;
+        }
 
-    return cache_.found(name);
+        return cache_.found(name);
+    }
+    else
+    {
+        return false;
+    }
 }
 
 
@@ -248,6 +256,7 @@ bool Foam::solution::read()
         if (dict.found("cache"))
         {
             cache_ = dict.subDict("cache");
+            caching_ = cache_.lookupOrDefault<Switch>("active", true);
         }
 
         if (dict.found("relaxationFactors"))
diff --git a/src/OpenFOAM/matrices/solution/solution.H b/src/OpenFOAM/matrices/solution/solution.H
index 36724bf495..ac1304fedc 100644
--- a/src/OpenFOAM/matrices/solution/solution.H
+++ b/src/OpenFOAM/matrices/solution/solution.H
@@ -37,6 +37,7 @@ SourceFiles
 #define solution_H
 
 #include "IOdictionary.H"
+#include "Switch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -56,6 +57,9 @@ class solution
         //- Dictionary of temporary fields to cache
         dictionary cache_;
 
+        //- Switch for the caching mechanism
+        Switch caching_;
+
         //- Dictionary of relaxation factors for all the fields
         dictionary relaxationFactors_;
 
diff --git a/src/OpenFOAM/meshes/primitiveShapes/objectHit/objectHit.H b/src/OpenFOAM/meshes/primitiveShapes/objectHit/objectHit.H
index fe2bf09865..16e1fc3a22 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/objectHit/objectHit.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/objectHit/objectHit.H
@@ -76,7 +76,7 @@ public:
         {}
 
         //- Construct from components
-        objectHit(const bool success, const label& obj)
+        objectHit(const bool success, const label obj)
         :
             hit_(success),
             hitObject_(obj)
@@ -111,7 +111,7 @@ public:
         {
             return ((a.hit_ == b.hit_) && (a.hitObject_ == b.hitObject_));
         }
-            
+
         friend bool operator!=(const objectHit& a, const objectHit& b)
         {
             return (!(a == b));
diff --git a/src/OpenFOAM/primitives/hashes/Hash/Hash.H b/src/OpenFOAM/primitives/hashes/Hash/Hash.H
index dd16517dfe..b998daea1b 100644
--- a/src/OpenFOAM/primitives/hashes/Hash/Hash.H
+++ b/src/OpenFOAM/primitives/hashes/Hash/Hash.H
@@ -107,7 +107,7 @@ public:
     //- Incrementally hash a label.
     //  This will necessarily return a different value than the
     //  non-incremental version.
-    unsigned operator()(const label& p, unsigned seed) const
+    unsigned operator()(const label p, unsigned seed) const
     {
         return Hasher(&p, sizeof(label), seed);
     }
@@ -115,11 +115,10 @@ public:
     //- Return the unsigned representation of a label.
     //  This helps if people have relied on the hash value corresponding to
     //  the natural order.
-    unsigned operator()(const label& p) const
+    unsigned operator()(const label p) const
     {
         return p;
     }
-
 };
 
 
diff --git a/src/OpenFOAM/primitives/random/Random.C b/src/OpenFOAM/primitives/random/Random.C
index 848b3a5547..570817e692 100644
--- a/src/OpenFOAM/primitives/random/Random.C
+++ b/src/OpenFOAM/primitives/random/Random.C
@@ -47,8 +47,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-// construct given seed
-Random::Random(const label& seed)
+Random::Random(const label seed)
 {
     if (seed > 1)
     {
diff --git a/src/OpenFOAM/primitives/random/Random.H b/src/OpenFOAM/primitives/random/Random.H
index f611f90600..2339e607b5 100644
--- a/src/OpenFOAM/primitives/random/Random.H
+++ b/src/OpenFOAM/primitives/random/Random.H
@@ -61,7 +61,7 @@ public:
     // Constructors
 
         //- Construct given seed
-        Random(const label&);
+        Random(const label);
 
 
     // Member functions
@@ -70,14 +70,19 @@ public:
 
         //- scalar [0..1] (so including 0,1)
         scalar scalar01();
+
         //- vector with every component scalar01
         vector vector01();
+
         //- sphericalTensor with every component scalar01
         sphericalTensor sphericalTensor01();
+
         //- symmTensor with every component scalar01
         symmTensor symmTensor01();
+
         //- tensor with every component scalar01
         tensor tensor01();
+
         //- label [lower..upper]
         label integer(const label lower, const label upper);
 
diff --git a/src/conversion/ensight/file/ensightFile.C b/src/conversion/ensight/file/ensightFile.C
index 8cb28f8c5e..095d2659f8 100644
--- a/src/conversion/ensight/file/ensightFile.C
+++ b/src/conversion/ensight/file/ensightFile.C
@@ -74,7 +74,7 @@ bool Foam::ensightFile::allowUndef(bool value)
 }
 
 
-Foam::scalar Foam::ensightFile::undefValue(const scalar& value)
+Foam::scalar Foam::ensightFile::undefValue(const scalar value)
 {
     // enable its use too
     allowUndef_ = true;
@@ -133,7 +133,7 @@ Foam::Ostream& Foam::ensightFile::write(const string& value)
 }
 
 
-Foam::Ostream& Foam::ensightFile::write(const label& value)
+Foam::Ostream& Foam::ensightFile::write(const label value)
 {
     if (format() == IOstream::BINARY)
     {
@@ -157,7 +157,7 @@ Foam::Ostream& Foam::ensightFile::write(const label& value)
 
 Foam::Ostream& Foam::ensightFile::write
 (
-    const label& value,
+    const label value,
     const label fieldWidth
 )
 {
@@ -181,7 +181,7 @@ Foam::Ostream& Foam::ensightFile::write
 }
 
 
-Foam::Ostream& Foam::ensightFile::write(const scalar& value)
+Foam::Ostream& Foam::ensightFile::write(const scalar value)
 {
     if (format() == IOstream::BINARY)
     {
diff --git a/src/conversion/ensight/file/ensightFile.H b/src/conversion/ensight/file/ensightFile.H
index 1fa91f8561..ca62114160 100644
--- a/src/conversion/ensight/file/ensightFile.H
+++ b/src/conversion/ensight/file/ensightFile.H
@@ -103,7 +103,7 @@ public:
         //- Assign the value to represent undef in the results
         //  Returns the previous value
         //  NB: do not use values larger than floatScalarVGREAT
-        static scalar undefValue(const scalar&);
+        static scalar undefValue(const scalar);
 
 
     // Output
@@ -124,13 +124,13 @@ public:
         Ostream& write(const string& value);
 
         //- write integer as "%10d" or as binary
-        Ostream& write(const label& value);
+        Ostream& write(const label value);
 
         //- write integer with specified width or as binary
-        Ostream& write(const label& value, const label fieldWidth);
+        Ostream& write(const label value, const label fieldWidth);
 
         //- write float as "%12.5e" or as binary
-        Ostream& write(const scalar& value);
+        Ostream& write(const scalar value);
 
         //- Add carriage return to ascii stream
         void newline();
diff --git a/src/conversion/meshTables/cellTable.C b/src/conversion/meshTables/cellTable.C
index 6369f98d9d..a07cbaf96d 100644
--- a/src/conversion/meshTables/cellTable.C
+++ b/src/conversion/meshTables/cellTable.C
@@ -22,8 +22,6 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
 #include "cellTable.H"
@@ -81,7 +79,7 @@ void Foam::cellTable::addDefaults()
 
 void Foam::cellTable::setEntry
 (
-    const label& id,
+    const label id,
     const word& keyWord,
     const word& value
 )
@@ -192,7 +190,7 @@ Foam::Map<Foam::word> Foam::cellTable::names
 }
 
 
-Foam::word Foam::cellTable::name(const label& id) const
+Foam::word Foam::cellTable::name(const label id) const
 {
     word theName("cellTable_" + Foam::name(id));
 
@@ -289,19 +287,19 @@ Foam::Map<Foam::word> Foam::cellTable::shells() const
 
 
 
-void Foam::cellTable::setMaterial(const label& id, const word& matlType)
+void Foam::cellTable::setMaterial(const label id, const word& matlType)
 {
     setEntry(id, "MaterialType", matlType);
 }
 
 
-void Foam::cellTable::setName(const label& id, const word& name)
+void Foam::cellTable::setName(const label id, const word& name)
 {
     setEntry(id, "Label", name);
 }
 
 
-void Foam::cellTable::setName(const label& id)
+void Foam::cellTable::setName(const label id)
 {
     iterator iter = find(id);
 
diff --git a/src/conversion/meshTables/cellTable.H b/src/conversion/meshTables/cellTable.H
index 4d0079f032..ff8c1fafb6 100644
--- a/src/conversion/meshTables/cellTable.H
+++ b/src/conversion/meshTables/cellTable.H
@@ -96,7 +96,7 @@ class cellTable
         //- Add required entries - MaterialType
         void addDefaults();
 
-        void setEntry(const label& id, const word& keyWord, const word& value);
+        void setEntry(const label id, const word& keyWord, const word& value);
 
         //- Disallow default bitwise copy construct
         cellTable(const cellTable&);
@@ -133,7 +133,7 @@ public:
 
         //- Return the name corresponding to id
         //  returns cellTable_ID if not otherwise defined
-        word name(const label& id) const;
+        word name(const label id) const;
 
         //- Return a Map of (id => name)
         Map<word> names() const;
@@ -157,13 +157,13 @@ public:
         Map<word> materialTypes() const;
 
         //- Assign material Type
-        void setMaterial(const label&, const word&);
+        void setMaterial(const label, const word&);
 
         //- Assign name
-        void setName(const label&, const word&);
+        void setName(const label, const word&);
 
         //- Assign default name if not already set
-        void setName(const label&);
+        void setName(const label);
 
         //- Read constant/cellTable
         void readDict
diff --git a/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C b/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C
index 7820bcd70b..2be5111187 100644
--- a/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C
+++ b/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C
@@ -1536,7 +1536,7 @@ void Foam::faceCoupleInfo::perfectPointMatch
         FatalErrorIn
         (
             "faceCoupleInfo::perfectPointMatch"
-            "(const scalar&, const bool)"
+            "(const scalar, const bool)"
         )   << "Did not match all of the master faces to the slave faces"
             << endl
             << "This usually means that the slave patch and master patch"
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C
index 926f2f4913..f338c6cf06 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C
@@ -55,7 +55,7 @@ namespace Foam
     class ifEqEqOp
     {
         public:
-        void operator()(label& x, const label& y) const
+        void operator()(label x, const label y) const
         {
             x = (x==y) ? x : value;
         }
diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index 7ecbf878a6..c387ac994d 100644
--- a/src/finiteVolume/Make/files
+++ b/src/finiteVolume/Make/files
@@ -285,6 +285,7 @@ $(divSchemes)/gaussDivScheme/gaussDivSchemes.C
 gradSchemes = finiteVolume/gradSchemes
 $(gradSchemes)/gradScheme/gradSchemes.C
 $(gradSchemes)/gaussGrad/gaussGrads.C
+
 $(gradSchemes)/leastSquaresGrad/leastSquaresVectors.C
 $(gradSchemes)/leastSquaresGrad/leastSquaresGrads.C
 $(gradSchemes)/extendedLeastSquaresGrad/extendedLeastSquaresVectors.C
diff --git a/src/finiteVolume/finiteVolume/convectionSchemes/convectionScheme/convectionScheme.H b/src/finiteVolume/finiteVolume/convectionSchemes/convectionScheme/convectionScheme.H
index c1e46a9e10..25dcbeb130 100644
--- a/src/finiteVolume/finiteVolume/convectionSchemes/convectionScheme/convectionScheme.H
+++ b/src/finiteVolume/finiteVolume/convectionSchemes/convectionScheme/convectionScheme.H
@@ -177,7 +177,7 @@ public:
         virtual tmp<fvMatrix<Type> > fvmDiv
         (
             const surfaceScalarField&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         ) const = 0;
 
         virtual tmp<GeometricField<Type, fvPatchField, volMesh> > fvcDiv
diff --git a/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.C b/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.C
index a219d80281..e2077d75ef 100644
--- a/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.C
+++ b/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.C
@@ -69,7 +69,7 @@ tmp<fvMatrix<Type> >
 gaussConvectionScheme<Type>::fvmDiv
 (
     const surfaceScalarField& faceFlux,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 ) const
 {
     tmp<surfaceScalarField> tweights = tinterpScheme_().weights(vf);
@@ -89,9 +89,9 @@ gaussConvectionScheme<Type>::fvmDiv
     fvm.upper() = fvm.lower() + faceFlux.internalField();
     fvm.negSumDiag();
 
-    forAll(fvm.psi().boundaryField(), patchI)
+    forAll(vf.boundaryField(), patchI)
     {
-        const fvPatchField<Type>& psf = fvm.psi().boundaryField()[patchI];
+        const fvPatchField<Type>& psf = vf.boundaryField()[patchI];
         const fvsPatchScalarField& patchFlux = faceFlux.boundaryField()[patchI];
         const fvsPatchScalarField& pw = weights.boundaryField()[patchI];
 
diff --git a/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.H b/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.H
index 7a4e5da155..8202f270dc 100644
--- a/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.H
+++ b/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.H
@@ -124,7 +124,7 @@ public:
         tmp<fvMatrix<Type> > fvmDiv
         (
             const surfaceScalarField&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         ) const;
 
         tmp<GeometricField<Type, fvPatchField, volMesh> > fvcDiv
diff --git a/src/finiteVolume/finiteVolume/convectionSchemes/multivariateGaussConvectionScheme/multivariateGaussConvectionScheme.C b/src/finiteVolume/finiteVolume/convectionSchemes/multivariateGaussConvectionScheme/multivariateGaussConvectionScheme.C
index d147d461ee..61cbba5ceb 100644
--- a/src/finiteVolume/finiteVolume/convectionSchemes/multivariateGaussConvectionScheme/multivariateGaussConvectionScheme.C
+++ b/src/finiteVolume/finiteVolume/convectionSchemes/multivariateGaussConvectionScheme/multivariateGaussConvectionScheme.C
@@ -22,8 +22,6 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-    
 \*---------------------------------------------------------------------------*/
 
 #include "multivariateGaussConvectionScheme.H"
@@ -81,7 +79,7 @@ tmp<fvMatrix<Type> >
 multivariateGaussConvectionScheme<Type>::fvmDiv
 (
     const surfaceScalarField& faceFlux,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 ) const
 {
     return gaussConvectionScheme<Type>
diff --git a/src/finiteVolume/finiteVolume/convectionSchemes/multivariateGaussConvectionScheme/multivariateGaussConvectionScheme.H b/src/finiteVolume/finiteVolume/convectionSchemes/multivariateGaussConvectionScheme/multivariateGaussConvectionScheme.H
index b8ccb25ed5..419628ddfa 100644
--- a/src/finiteVolume/finiteVolume/convectionSchemes/multivariateGaussConvectionScheme/multivariateGaussConvectionScheme.H
+++ b/src/finiteVolume/finiteVolume/convectionSchemes/multivariateGaussConvectionScheme/multivariateGaussConvectionScheme.H
@@ -114,7 +114,7 @@ public:
         tmp<fvMatrix<Type> > fvmDiv
         (
             const surfaceScalarField&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         ) const;
 
         tmp<GeometricField<Type, fvPatchField, volMesh> > fvcDiv
diff --git a/src/finiteVolume/finiteVolume/d2dt2Schemes/EulerD2dt2Scheme/EulerD2dt2Scheme.C b/src/finiteVolume/finiteVolume/d2dt2Schemes/EulerD2dt2Scheme/EulerD2dt2Scheme.C
index abc203a73b..0285084dff 100644
--- a/src/finiteVolume/finiteVolume/d2dt2Schemes/EulerD2dt2Scheme/EulerD2dt2Scheme.C
+++ b/src/finiteVolume/finiteVolume/d2dt2Schemes/EulerD2dt2Scheme/EulerD2dt2Scheme.C
@@ -21,7 +21,7 @@ License
     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 "EulerD2dt2Scheme.H"
@@ -181,7 +181,7 @@ EulerD2dt2Scheme<Type>::fvcD2dt2
                     coefft
                    *(rho.boundaryField() + rho.oldTime().boundaryField())
                    *vf.boundaryField()
-                        
+
                   - (
                         coefft
                        *(
@@ -232,7 +232,7 @@ template<class Type>
 tmp<fvMatrix<Type> >
 EulerD2dt2Scheme<Type>::fvmD2dt2
 (
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > tfvm
@@ -292,7 +292,7 @@ tmp<fvMatrix<Type> >
 EulerD2dt2Scheme<Type>::fvmD2dt2
 (
     const dimensionedScalar& rho,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > tfvm
@@ -353,7 +353,7 @@ tmp<fvMatrix<Type> >
 EulerD2dt2Scheme<Type>::fvmD2dt2
 (
     const volScalarField& rho,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > tfvm
diff --git a/src/finiteVolume/finiteVolume/d2dt2Schemes/EulerD2dt2Scheme/EulerD2dt2Scheme.H b/src/finiteVolume/finiteVolume/d2dt2Schemes/EulerD2dt2Scheme/EulerD2dt2Scheme.H
index bca75562ed..0823f225dc 100644
--- a/src/finiteVolume/finiteVolume/d2dt2Schemes/EulerD2dt2Scheme/EulerD2dt2Scheme.H
+++ b/src/finiteVolume/finiteVolume/d2dt2Schemes/EulerD2dt2Scheme/EulerD2dt2Scheme.H
@@ -109,19 +109,19 @@ public:
 
         tmp<fvMatrix<Type> > fvmD2dt2
         (
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         tmp<fvMatrix<Type> > fvmD2dt2
         (
             const dimensionedScalar&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         tmp<fvMatrix<Type> > fvmD2dt2
         (
             const volScalarField&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 };
 
diff --git a/src/finiteVolume/finiteVolume/d2dt2Schemes/d2dt2Scheme/d2dt2Scheme.H b/src/finiteVolume/finiteVolume/d2dt2Schemes/d2dt2Scheme/d2dt2Scheme.H
index 9e5668ecf3..452373b501 100644
--- a/src/finiteVolume/finiteVolume/d2dt2Schemes/d2dt2Scheme/d2dt2Scheme.H
+++ b/src/finiteVolume/finiteVolume/d2dt2Schemes/d2dt2Scheme/d2dt2Scheme.H
@@ -153,19 +153,19 @@ public:
 
         virtual tmp<fvMatrix<Type> > fvmD2dt2
         (
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         ) = 0;
 
         virtual tmp<fvMatrix<Type> > fvmD2dt2
         (
             const dimensionedScalar&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         ) = 0;
 
         virtual tmp<fvMatrix<Type> > fvmD2dt2
         (
             const volScalarField&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         ) = 0;
 };
 
diff --git a/src/finiteVolume/finiteVolume/d2dt2Schemes/steadyStateD2dt2Scheme/steadyStateD2dt2Scheme.C b/src/finiteVolume/finiteVolume/d2dt2Schemes/steadyStateD2dt2Scheme/steadyStateD2dt2Scheme.C
index 121df11087..9d810c401c 100644
--- a/src/finiteVolume/finiteVolume/d2dt2Schemes/steadyStateD2dt2Scheme/steadyStateD2dt2Scheme.C
+++ b/src/finiteVolume/finiteVolume/d2dt2Schemes/steadyStateD2dt2Scheme/steadyStateD2dt2Scheme.C
@@ -21,7 +21,7 @@ License
     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 "steadyStateD2dt2Scheme.H"
@@ -107,7 +107,7 @@ template<class Type>
 tmp<fvMatrix<Type> >
 steadyStateD2dt2Scheme<Type>::fvmD2dt2
 (
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > tfvm
@@ -128,7 +128,7 @@ tmp<fvMatrix<Type> >
 steadyStateD2dt2Scheme<Type>::fvmD2dt2
 (
     const dimensionedScalar& rho,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > tfvm
@@ -149,7 +149,7 @@ tmp<fvMatrix<Type> >
 steadyStateD2dt2Scheme<Type>::fvmD2dt2
 (
     const volScalarField& rho,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > tfvm
diff --git a/src/finiteVolume/finiteVolume/d2dt2Schemes/steadyStateD2dt2Scheme/steadyStateD2dt2Scheme.H b/src/finiteVolume/finiteVolume/d2dt2Schemes/steadyStateD2dt2Scheme/steadyStateD2dt2Scheme.H
index 004fb22a5e..9b9327ceee 100644
--- a/src/finiteVolume/finiteVolume/d2dt2Schemes/steadyStateD2dt2Scheme/steadyStateD2dt2Scheme.H
+++ b/src/finiteVolume/finiteVolume/d2dt2Schemes/steadyStateD2dt2Scheme/steadyStateD2dt2Scheme.H
@@ -108,19 +108,19 @@ public:
 
         tmp<fvMatrix<Type> > fvmD2dt2
         (
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         tmp<fvMatrix<Type> > fvmD2dt2
         (
             const dimensionedScalar&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         tmp<fvMatrix<Type> > fvmD2dt2
         (
             const volScalarField&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 };
 
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C
index 6928e20a2e..a63a644427 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C
@@ -21,7 +21,7 @@ License
     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 "CoEulerDdtScheme.H"
@@ -69,10 +69,10 @@ tmp<volScalarField> CoEulerDdtScheme<Type>::CorDeltaT() const
 
     forAll(owner, faceI)
     {
-        corDeltaT[owner[faceI]] = 
+        corDeltaT[owner[faceI]] =
             max(corDeltaT[owner[faceI]], cofrDeltaT[faceI]);
 
-        corDeltaT[neighbour[faceI]] = 
+        corDeltaT[neighbour[faceI]] =
             max(corDeltaT[neighbour[faceI]], cofrDeltaT[faceI]);
     }
 
@@ -127,7 +127,7 @@ tmp<surfaceScalarField> CoEulerDdtScheme<Type>::CofrDeltaT() const
         const volScalarField& rho =
             static_cast<const objectRegistry&>(mesh())
            .lookupObject<volScalarField>(rhoName_).oldTime();
- 
+
         surfaceScalarField Co
         (
             mesh().surfaceInterpolation::deltaCoeffs()
@@ -369,7 +369,7 @@ template<class Type>
 tmp<fvMatrix<Type> >
 CoEulerDdtScheme<Type>::fvmDdt
 (
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > tfvm
@@ -386,7 +386,7 @@ CoEulerDdtScheme<Type>::fvmDdt
     scalarField rDeltaT = CorDeltaT()().internalField();
 
     fvm.diag() = rDeltaT*mesh().V();
-    
+
     if (mesh().moving())
     {
         fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().V0();
@@ -405,7 +405,7 @@ tmp<fvMatrix<Type> >
 CoEulerDdtScheme<Type>::fvmDdt
 (
     const dimensionedScalar& rho,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > tfvm
@@ -421,7 +421,7 @@ CoEulerDdtScheme<Type>::fvmDdt
     scalarField rDeltaT = CorDeltaT()().internalField();
 
     fvm.diag() = rDeltaT*rho.value()*mesh().V();
-    
+
     if (mesh().moving())
     {
         fvm.source() = rDeltaT
@@ -442,7 +442,7 @@ tmp<fvMatrix<Type> >
 CoEulerDdtScheme<Type>::fvmDdt
 (
     const volScalarField& rho,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > tfvm
@@ -588,7 +588,7 @@ CoEulerDdtScheme<Type>::fvcDdtPhiCorr
                 )
             );
         }
-        else if 
+        else if
         (
             U.dimensions() == dimVelocity
          && phi.dimensions() == dimDensity*dimVelocity*dimArea
@@ -617,7 +617,7 @@ CoEulerDdtScheme<Type>::fvcDdtPhiCorr
                 )
             );
         }
-        else if 
+        else if
         (
             U.dimensions() == dimDensity*dimVelocity
          && phi.dimensions() == dimDensity*dimVelocity*dimArea
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.H
index f2b54169ed..65cd320d2d 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.H
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.H
@@ -142,19 +142,19 @@ public:
 
         tmp<fvMatrix<Type> > fvmDdt
         (
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         tmp<fvMatrix<Type> > fvmDdt
         (
             const dimensionedScalar&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         tmp<fvMatrix<Type> > fvmDdt
         (
             const volScalarField&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         typedef typename ddtScheme<Type>::fluxFieldType fluxFieldType;
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicholsonDdtScheme/CrankNicholsonDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicholsonDdtScheme/CrankNicholsonDdtScheme.C
index 027f0a2b30..172bea60aa 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicholsonDdtScheme/CrankNicholsonDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicholsonDdtScheme/CrankNicholsonDdtScheme.C
@@ -625,7 +625,7 @@ template<class Type>
 tmp<fvMatrix<Type> >
 CrankNicholsonDdtScheme<Type>::fvmDdt
 (
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     DDt0Field<GeometricField<Type, fvPatchField, volMesh> >& ddt0 =
@@ -709,7 +709,7 @@ tmp<fvMatrix<Type> >
 CrankNicholsonDdtScheme<Type>::fvmDdt
 (
     const dimensionedScalar& rho,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     DDt0Field<GeometricField<Type, fvPatchField, volMesh> >& ddt0 =
@@ -791,7 +791,7 @@ tmp<fvMatrix<Type> >
 CrankNicholsonDdtScheme<Type>::fvmDdt
 (
     const volScalarField& rho,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     DDt0Field<GeometricField<Type, fvPatchField, volMesh> >& ddt0 =
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicholsonDdtScheme/CrankNicholsonDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicholsonDdtScheme/CrankNicholsonDdtScheme.H
index 87ca9a086d..d3177f6a95 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicholsonDdtScheme/CrankNicholsonDdtScheme.H
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicholsonDdtScheme/CrankNicholsonDdtScheme.H
@@ -90,7 +90,7 @@ class CrankNicholsonDdtScheme
 
             //- Return the start-time index
             label startTimeIndex() const;
-            
+
             //- Cast to the underlying GeoField
             GeoField& operator()();
 
@@ -213,19 +213,19 @@ public:
 
         tmp<fvMatrix<Type> > fvmDdt
         (
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         tmp<fvMatrix<Type> > fvmDdt
         (
             const dimensionedScalar&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         tmp<fvMatrix<Type> > fvmDdt
         (
             const volScalarField&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         typedef typename ddtScheme<Type>::fluxFieldType fluxFieldType;
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C
index 5b43623c61..91e4cb403a 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C
@@ -262,7 +262,7 @@ template<class Type>
 tmp<fvMatrix<Type> >
 EulerDdtScheme<Type>::fvmDdt
 (
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > tfvm
@@ -298,7 +298,7 @@ tmp<fvMatrix<Type> >
 EulerDdtScheme<Type>::fvmDdt
 (
     const dimensionedScalar& rho,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > tfvm
@@ -335,7 +335,7 @@ tmp<fvMatrix<Type> >
 EulerDdtScheme<Type>::fvmDdt
 (
     const volScalarField& rho,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > tfvm
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.H
index b06ef32761..3d38050955 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.H
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.H
@@ -120,19 +120,19 @@ public:
 
         tmp<fvMatrix<Type> > fvmDdt
         (
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         tmp<fvMatrix<Type> > fvmDdt
         (
             const dimensionedScalar&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         tmp<fvMatrix<Type> > fvmDdt
         (
             const volScalarField&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         typedef typename ddtScheme<Type>::fluxFieldType fluxFieldType;
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C
index a06e30e725..4906f63567 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C
@@ -21,7 +21,7 @@ License
     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 "SLTSDdtScheme.H"
@@ -369,7 +369,7 @@ template<class Type>
 tmp<fvMatrix<Type> >
 SLTSDdtScheme<Type>::fvmDdt
 (
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > tfvm
@@ -388,7 +388,7 @@ SLTSDdtScheme<Type>::fvmDdt
     Info<< "max/min rDeltaT " << max(rDeltaT) << " " << min(rDeltaT) << endl;
 
     fvm.diag() = rDeltaT*mesh().V();
-    
+
     if (mesh().moving())
     {
         fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().V0();
@@ -407,7 +407,7 @@ tmp<fvMatrix<Type> >
 SLTSDdtScheme<Type>::fvmDdt
 (
     const dimensionedScalar& rho,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > tfvm
@@ -423,7 +423,7 @@ SLTSDdtScheme<Type>::fvmDdt
     scalarField rDeltaT = SLrDeltaT()().internalField();
 
     fvm.diag() = rDeltaT*rho.value()*mesh().V();
-    
+
     if (mesh().moving())
     {
         fvm.source() = rDeltaT
@@ -444,7 +444,7 @@ tmp<fvMatrix<Type> >
 SLTSDdtScheme<Type>::fvmDdt
 (
     const volScalarField& rho,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > tfvm
@@ -590,7 +590,7 @@ SLTSDdtScheme<Type>::fvcDdtPhiCorr
                 )
             );
         }
-        else if 
+        else if
         (
             U.dimensions() == dimVelocity
          && phi.dimensions() == dimDensity*dimVelocity*dimArea
@@ -619,7 +619,7 @@ SLTSDdtScheme<Type>::fvcDdtPhiCorr
                 )
             );
         }
-        else if 
+        else if
         (
             U.dimensions() == dimDensity*dimVelocity
          && phi.dimensions() == dimDensity*dimVelocity*dimArea
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.H
index 850f358299..f300297e6c 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.H
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.H
@@ -143,19 +143,19 @@ public:
 
         tmp<fvMatrix<Type> > fvmDdt
         (
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         tmp<fvMatrix<Type> > fvmDdt
         (
             const dimensionedScalar&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         tmp<fvMatrix<Type> > fvmDdt
         (
             const volScalarField&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         typedef typename ddtScheme<Type>::fluxFieldType fluxFieldType;
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C
index 5b55f27a1e..6051903a9c 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C
@@ -361,7 +361,7 @@ template<class Type>
 tmp<fvMatrix<Type> >
 backwardDdtScheme<Type>::fvmDdt
 (
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > tfvm
@@ -413,7 +413,7 @@ tmp<fvMatrix<Type> >
 backwardDdtScheme<Type>::fvmDdt
 (
     const dimensionedScalar& rho,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > tfvm
@@ -464,7 +464,7 @@ tmp<fvMatrix<Type> >
 backwardDdtScheme<Type>::fvmDdt
 (
     const volScalarField& rho,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > tfvm
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.H
index 374d05031a..7aceb79c89 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.H
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.H
@@ -131,19 +131,19 @@ public:
 
         tmp<fvMatrix<Type> > fvmDdt
         (
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         tmp<fvMatrix<Type> > fvmDdt
         (
             const dimensionedScalar&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         tmp<fvMatrix<Type> > fvmDdt
         (
             const volScalarField&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         typedef typename ddtScheme<Type>::fluxFieldType fluxFieldType;
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/boundedBackwardDdtScheme/boundedBackwardDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/boundedBackwardDdtScheme/boundedBackwardDdtScheme.C
index c7b8d49866..236f66cc29 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/boundedBackwardDdtScheme/boundedBackwardDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/boundedBackwardDdtScheme/boundedBackwardDdtScheme.C
@@ -21,7 +21,7 @@ License
     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 "boundedBackwardDdtScheme.H"
@@ -413,7 +413,7 @@ boundedBackwardDdtScheme::fvcDdt
 tmp<fvScalarMatrix>
 boundedBackwardDdtScheme::fvmDdt
 (
-    volScalarField& vf
+    const volScalarField& vf
 )
 {
     tmp<fvScalarMatrix> tfvm
@@ -484,7 +484,7 @@ tmp<fvScalarMatrix>
 boundedBackwardDdtScheme::fvmDdt
 (
     const dimensionedScalar& rho,
-    volScalarField& vf
+    const volScalarField& vf
 )
 {
     tmp<fvScalarMatrix> tfvm
@@ -554,7 +554,7 @@ tmp<fvScalarMatrix>
 boundedBackwardDdtScheme::fvmDdt
 (
     const volScalarField& rho,
-    volScalarField& vf
+    const volScalarField& vf
 )
 {
     tmp<fvScalarMatrix> tfvm
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/boundedBackwardDdtScheme/boundedBackwardDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/boundedBackwardDdtScheme/boundedBackwardDdtScheme.H
index 7bff042a75..f1b4121f6e 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/boundedBackwardDdtScheme/boundedBackwardDdtScheme.H
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/boundedBackwardDdtScheme/boundedBackwardDdtScheme.H
@@ -142,19 +142,19 @@ public:
 
         tmp<fvScalarMatrix> fvmDdt
         (
-            volScalarField&
+            const volScalarField&
         );
 
         tmp<fvScalarMatrix> fvmDdt
         (
             const dimensionedScalar&,
-            volScalarField&
+            const volScalarField&
         );
 
         tmp<fvScalarMatrix> fvmDdt
         (
             const volScalarField&,
-            volScalarField&
+            const volScalarField&
         );
 
         tmp<surfaceScalarField> fvcDdtPhiCorr
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/boundedBackwardDdtScheme/boundedBackwardDdtSchemes.C b/src/finiteVolume/finiteVolume/ddtSchemes/boundedBackwardDdtScheme/boundedBackwardDdtSchemes.C
index 666cb2be3b..67a1d1f69c 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/boundedBackwardDdtScheme/boundedBackwardDdtSchemes.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/boundedBackwardDdtScheme/boundedBackwardDdtSchemes.C
@@ -21,7 +21,7 @@ License
     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 "boundedBackwardDdtScheme.H"
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.H
index 31f8fa090d..e454e9409a 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.H
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.H
@@ -164,19 +164,19 @@ public:
 
         virtual tmp<fvMatrix<Type> > fvmDdt
         (
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         ) = 0;
 
         virtual tmp<fvMatrix<Type> > fvmDdt
         (
             const dimensionedScalar&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         ) = 0;
 
         virtual tmp<fvMatrix<Type> > fvmDdt
         (
             const volScalarField&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         ) = 0;
 
 
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.C
index 7bd7b6f8a8..3b5736adb8 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.C
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.C
@@ -21,7 +21,7 @@ License
     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 "steadyStateDdtScheme.H"
@@ -162,7 +162,7 @@ template<class Type>
 tmp<fvMatrix<Type> >
 steadyStateDdtScheme<Type>::fvmDdt
 (
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > tfvm
@@ -183,7 +183,7 @@ tmp<fvMatrix<Type> >
 steadyStateDdtScheme<Type>::fvmDdt
 (
     const dimensionedScalar& rho,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > tfvm
@@ -204,7 +204,7 @@ tmp<fvMatrix<Type> >
 steadyStateDdtScheme<Type>::fvmDdt
 (
     const volScalarField& rho,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > tfvm
diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.H b/src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.H
index 54208441b4..0af9bb66db 100644
--- a/src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.H
+++ b/src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.H
@@ -119,19 +119,19 @@ public:
 
         tmp<fvMatrix<Type> > fvmDdt
         (
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         tmp<fvMatrix<Type> > fvmDdt
         (
             const dimensionedScalar&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         tmp<fvMatrix<Type> > fvmDdt
         (
             const volScalarField&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         typedef typename ddtScheme<Type>::fluxFieldType fluxFieldType;
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcGrad.C b/src/finiteVolume/finiteVolume/fvc/fvcGrad.C
index 9d61a45e40..9f5c2bd98f 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcGrad.C
+++ b/src/finiteVolume/finiteVolume/fvc/fvcGrad.C
@@ -54,7 +54,7 @@ grad
     const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf
 )
 {
-    return fv::gaussGrad<Type>::grad(ssf);
+    return fv::gaussGrad<Type>::gradf(ssf, "grad(" + ssf.name() + ')');
 }
 
 
@@ -99,7 +99,7 @@ grad
     (
         vf.mesh(),
         vf.mesh().gradScheme(name)
-    )().grad(vf);
+    )().grad(vf, name);
 }
 
 
diff --git a/src/finiteVolume/finiteVolume/fvm/fvmD2dt2.C b/src/finiteVolume/finiteVolume/fvm/fvmD2dt2.C
index 4df081cd07..81f37cd099 100644
--- a/src/finiteVolume/finiteVolume/fvm/fvmD2dt2.C
+++ b/src/finiteVolume/finiteVolume/fvm/fvmD2dt2.C
@@ -22,9 +22,6 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-    
-
 \*---------------------------------------------------------------------------*/
 
 #include "volFields.H"
@@ -48,7 +45,7 @@ template<class Type>
 tmp<fvMatrix<Type> >
 d2dt2
 (
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     return fv::d2dt2Scheme<Type>::New
@@ -64,7 +61,7 @@ tmp<fvMatrix<Type> >
 d2dt2
 (
     const dimensionedScalar& rho,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     return fv::d2dt2Scheme<Type>::New
@@ -80,7 +77,7 @@ tmp<fvMatrix<Type> >
 d2dt2
 (
     const volScalarField& rho,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     return fv::d2dt2Scheme<Type>::New
diff --git a/src/finiteVolume/finiteVolume/fvm/fvmD2dt2.H b/src/finiteVolume/finiteVolume/fvm/fvmD2dt2.H
index b322391ade..3d7e2a716b 100644
--- a/src/finiteVolume/finiteVolume/fvm/fvmD2dt2.H
+++ b/src/finiteVolume/finiteVolume/fvm/fvmD2dt2.H
@@ -54,20 +54,20 @@ namespace fvm
     tmp<fvMatrix<Type> > d2dt2
     (
         const dimensionedScalar&,
-        GeometricField<Type, fvPatchField, volMesh>&
+        const GeometricField<Type, fvPatchField, volMesh>&
     );
 
     template<class Type>
     tmp<fvMatrix<Type> > d2dt2
     (
-        GeometricField<Type, fvPatchField, volMesh>&
+        const GeometricField<Type, fvPatchField, volMesh>&
     );
 
     template<class Type>
     tmp<fvMatrix<Type> > d2dt2
     (
         const volScalarField&,
-        GeometricField<Type, fvPatchField, volMesh>&
+        const GeometricField<Type, fvPatchField, volMesh>&
     );
 }
 
diff --git a/src/finiteVolume/finiteVolume/fvm/fvmDdt.C b/src/finiteVolume/finiteVolume/fvm/fvmDdt.C
index c2e3e976af..d7610206d0 100644
--- a/src/finiteVolume/finiteVolume/fvm/fvmDdt.C
+++ b/src/finiteVolume/finiteVolume/fvm/fvmDdt.C
@@ -45,7 +45,7 @@ template<class Type>
 tmp<fvMatrix<Type> >
 ddt
 (
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     return fv::ddtScheme<Type>::New
@@ -61,7 +61,7 @@ tmp<fvMatrix<Type> >
 ddt
 (
     const oneField&,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     return ddt(vf);
@@ -73,7 +73,7 @@ tmp<fvMatrix<Type> >
 ddt
 (
     const dimensionedScalar& rho,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     return fv::ddtScheme<Type>::New
@@ -89,7 +89,7 @@ tmp<fvMatrix<Type> >
 ddt
 (
     const volScalarField& rho,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     return fv::ddtScheme<Type>::New
diff --git a/src/finiteVolume/finiteVolume/fvm/fvmDdt.H b/src/finiteVolume/finiteVolume/fvm/fvmDdt.H
index 5472f91f37..2b8124e559 100644
--- a/src/finiteVolume/finiteVolume/fvm/fvmDdt.H
+++ b/src/finiteVolume/finiteVolume/fvm/fvmDdt.H
@@ -54,28 +54,28 @@ namespace fvm
     template<class Type>
     tmp<fvMatrix<Type> > ddt
     (
-        GeometricField<Type, fvPatchField, volMesh>&
+        const GeometricField<Type, fvPatchField, volMesh>&
     );
 
     template<class Type>
     tmp<fvMatrix<Type> > ddt
     (
         const oneField&,
-        GeometricField<Type, fvPatchField, volMesh>&
+        const GeometricField<Type, fvPatchField, volMesh>&
     );
 
     template<class Type>
     tmp<fvMatrix<Type> > ddt
     (
         const dimensionedScalar&,
-        GeometricField<Type, fvPatchField, volMesh>&
+        const GeometricField<Type, fvPatchField, volMesh>&
     );
 
     template<class Type>
     tmp<fvMatrix<Type> > ddt
     (
         const volScalarField&,
-        GeometricField<Type, fvPatchField, volMesh>&
+        const GeometricField<Type, fvPatchField, volMesh>&
     );
 }
 
diff --git a/src/finiteVolume/finiteVolume/fvm/fvmDiv.C b/src/finiteVolume/finiteVolume/fvm/fvmDiv.C
index 25af21ef2a..2cc8d0e31f 100644
--- a/src/finiteVolume/finiteVolume/fvm/fvmDiv.C
+++ b/src/finiteVolume/finiteVolume/fvm/fvmDiv.C
@@ -22,9 +22,6 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-    
-
 \*---------------------------------------------------------------------------*/
 
 #include "fvmDiv.H"
@@ -49,7 +46,7 @@ tmp<fvMatrix<Type> >
 div
 (
     const surfaceScalarField& flux,
-    GeometricField<Type, fvPatchField, volMesh>& vf,
+    const GeometricField<Type, fvPatchField, volMesh>& vf,
     const word& name
 )
 {
@@ -66,7 +63,7 @@ tmp<fvMatrix<Type> >
 div
 (
     const tmp<surfaceScalarField>& tflux,
-    GeometricField<Type, fvPatchField, volMesh>& vf,
+    const GeometricField<Type, fvPatchField, volMesh>& vf,
     const word& name
 )
 {
@@ -81,7 +78,7 @@ tmp<fvMatrix<Type> >
 div
 (
     const surfaceScalarField& flux,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     return fvm::div(flux, vf, "div("+flux.name()+','+vf.name()+')');
@@ -92,7 +89,7 @@ tmp<fvMatrix<Type> >
 div
 (
     const tmp<surfaceScalarField>& tflux,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > Div(fvm::div(tflux(), vf));
diff --git a/src/finiteVolume/finiteVolume/fvm/fvmDiv.H b/src/finiteVolume/finiteVolume/fvm/fvmDiv.H
index c15fc72dec..9e45a883da 100644
--- a/src/finiteVolume/finiteVolume/fvm/fvmDiv.H
+++ b/src/finiteVolume/finiteVolume/fvm/fvmDiv.H
@@ -56,7 +56,7 @@ namespace fvm
     tmp<fvMatrix<Type> > div
     (
         const surfaceScalarField&,
-        GeometricField<Type, fvPatchField, volMesh>&,
+        const GeometricField<Type, fvPatchField, volMesh>&,
         const word& name
     );
 
@@ -64,7 +64,7 @@ namespace fvm
     tmp<fvMatrix<Type> > div
     (
         const tmp<surfaceScalarField>&,
-        GeometricField<Type, fvPatchField, volMesh>&,
+        const GeometricField<Type, fvPatchField, volMesh>&,
         const word& name
     );
 
@@ -73,14 +73,14 @@ namespace fvm
     tmp<fvMatrix<Type> > div
     (
         const surfaceScalarField&,
-        GeometricField<Type, fvPatchField, volMesh>&
+        const GeometricField<Type, fvPatchField, volMesh>&
     );
 
     template<class Type>
     tmp<fvMatrix<Type> > div
     (
         const tmp<surfaceScalarField>&,
-        GeometricField<Type, fvPatchField, volMesh>&
+        const GeometricField<Type, fvPatchField, volMesh>&
     );
 }
 
diff --git a/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.C b/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.C
index 4a60eadd3d..501456da3c 100644
--- a/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.C
+++ b/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.C
@@ -45,7 +45,7 @@ template<class Type>
 tmp<fvMatrix<Type> >
 laplacian
 (
-    GeometricField<Type, fvPatchField, volMesh>& vf,
+    const GeometricField<Type, fvPatchField, volMesh>& vf,
     const word& name
 )
 {
@@ -70,7 +70,7 @@ template<class Type>
 tmp<fvMatrix<Type> >
 laplacian
 (
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     surfaceScalarField Gamma
@@ -100,7 +100,7 @@ tmp<fvMatrix<Type> >
 laplacian
 (
     const zeroField&,
-    GeometricField<Type, fvPatchField, volMesh>& vf,
+    const GeometricField<Type, fvPatchField, volMesh>& vf,
     const word& name
 )
 {
@@ -116,7 +116,7 @@ tmp<fvMatrix<Type> >
 laplacian
 (
     const zeroField&,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     return tmp<fvMatrix<Type> >
@@ -131,7 +131,7 @@ tmp<fvMatrix<Type> >
 laplacian
 (
     const oneField&,
-    GeometricField<Type, fvPatchField, volMesh>& vf,
+    const GeometricField<Type, fvPatchField, volMesh>& vf,
     const word& name
 )
 {
@@ -144,7 +144,7 @@ tmp<fvMatrix<Type> >
 laplacian
 (
     const oneField&,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     return fvm::laplacian(vf);
@@ -156,11 +156,11 @@ tmp<fvMatrix<Type> >
 laplacian
 (
     const dimensioned<GType>& gamma,
-    GeometricField<Type, fvPatchField, volMesh>& vf,
+    const GeometricField<Type, fvPatchField, volMesh>& vf,
     const word& name
 )
 {
-    GeometricField<GType, fvsPatchField, surfaceMesh> Gamma
+    const GeometricField<GType, fvsPatchField, surfaceMesh> Gamma
     (
         IOobject
         (
@@ -182,10 +182,10 @@ tmp<fvMatrix<Type> >
 laplacian
 (
     const dimensioned<GType>& gamma,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
-    GeometricField<GType, fvsPatchField, surfaceMesh> Gamma
+    const GeometricField<GType, fvsPatchField, surfaceMesh> Gamma
     (
         IOobject
         (
@@ -209,7 +209,7 @@ tmp<fvMatrix<Type> >
 laplacian
 (
     const GeometricField<GType, fvPatchField, volMesh>& gamma,
-    GeometricField<Type, fvPatchField, volMesh>& vf,
+    const GeometricField<Type, fvPatchField, volMesh>& vf,
     const word& name
 )
 {
@@ -226,7 +226,7 @@ tmp<fvMatrix<Type> >
 laplacian
 (
     const tmp<GeometricField<GType, fvPatchField, volMesh> >& tgamma,
-    GeometricField<Type, fvPatchField, volMesh>& vf,
+    const GeometricField<Type, fvPatchField, volMesh>& vf,
     const word& name
 )
 {
@@ -241,7 +241,7 @@ tmp<fvMatrix<Type> >
 laplacian
 (
     const GeometricField<GType, fvPatchField, volMesh>& gamma,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     return fvm::laplacian
@@ -258,7 +258,7 @@ tmp<fvMatrix<Type> >
 laplacian
 (
     const tmp<GeometricField<GType, fvPatchField, volMesh> >& tgamma,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > Laplacian(fvm::laplacian(tgamma(), vf));
@@ -274,7 +274,7 @@ tmp<fvMatrix<Type> >
 laplacian
 (
     const GeometricField<GType, fvsPatchField, surfaceMesh>& gamma,
-    GeometricField<Type, fvPatchField, volMesh>& vf,
+    const GeometricField<Type, fvPatchField, volMesh>& vf,
     const word& name
 )
 {
@@ -291,7 +291,7 @@ tmp<fvMatrix<Type> >
 laplacian
 (
     const tmp<GeometricField<GType, fvsPatchField, surfaceMesh> >& tgamma,
-    GeometricField<Type, fvPatchField, volMesh>& vf,
+    const GeometricField<Type, fvPatchField, volMesh>& vf,
     const word& name
 )
 {
@@ -306,7 +306,7 @@ tmp<fvMatrix<Type> >
 laplacian
 (
     const GeometricField<GType, fvsPatchField, surfaceMesh>& gamma,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     return fvm::laplacian
@@ -323,7 +323,7 @@ tmp<fvMatrix<Type> >
 laplacian
 (
     const tmp<GeometricField<GType, fvsPatchField, surfaceMesh> >& tGamma,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > tfvm(fvm::laplacian(tGamma(), vf));
diff --git a/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.H b/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.H
index cc06ec5008..c41774361f 100644
--- a/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.H
+++ b/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.H
@@ -55,14 +55,14 @@ namespace fvm
     template<class Type>
     tmp<fvMatrix<Type> > laplacian
     (
-        GeometricField<Type, fvPatchField, volMesh>&,
+        const GeometricField<Type, fvPatchField, volMesh>&,
         const word&
     );
 
     template<class Type>
     tmp<fvMatrix<Type> > laplacian
     (
-        GeometricField<Type, fvPatchField, volMesh>&
+        const GeometricField<Type, fvPatchField, volMesh>&
     );
 
 
@@ -70,7 +70,7 @@ namespace fvm
     tmp<fvMatrix<Type> > laplacian
     (
         const zeroField&,
-        GeometricField<Type, fvPatchField, volMesh>&,
+        const GeometricField<Type, fvPatchField, volMesh>&,
         const word&
     );
 
@@ -78,7 +78,7 @@ namespace fvm
     tmp<fvMatrix<Type> > laplacian
     (
         const zeroField&,
-        GeometricField<Type, fvPatchField, volMesh>&
+        const GeometricField<Type, fvPatchField, volMesh>&
     );
 
 
@@ -86,7 +86,7 @@ namespace fvm
     tmp<fvMatrix<Type> > laplacian
     (
         const oneField&,
-        GeometricField<Type, fvPatchField, volMesh>&,
+        const GeometricField<Type, fvPatchField, volMesh>&,
         const word&
     );
 
@@ -94,7 +94,7 @@ namespace fvm
     tmp<fvMatrix<Type> > laplacian
     (
         const oneField&,
-        GeometricField<Type, fvPatchField, volMesh>&
+        const GeometricField<Type, fvPatchField, volMesh>&
     );
 
 
@@ -102,7 +102,7 @@ namespace fvm
     tmp<fvMatrix<Type> > laplacian
     (
         const dimensioned<GType>&,
-        GeometricField<Type, fvPatchField, volMesh>&,
+        const GeometricField<Type, fvPatchField, volMesh>&,
         const word&
     );
 
@@ -110,7 +110,7 @@ namespace fvm
     tmp<fvMatrix<Type> > laplacian
     (
         const dimensioned<GType>&,
-        GeometricField<Type, fvPatchField, volMesh>&
+        const GeometricField<Type, fvPatchField, volMesh>&
     );
 
 
@@ -118,7 +118,7 @@ namespace fvm
     tmp<fvMatrix<Type> > laplacian
     (
         const GeometricField<GType, fvPatchField, volMesh>&,
-        GeometricField<Type, fvPatchField, volMesh>&,
+        const GeometricField<Type, fvPatchField, volMesh>&,
         const word&
     );
 
@@ -126,7 +126,7 @@ namespace fvm
     tmp<fvMatrix<Type> > laplacian
     (
         const GeometricField<GType, fvPatchField, volMesh>&,
-        GeometricField<Type, fvPatchField, volMesh>&
+        const GeometricField<Type, fvPatchField, volMesh>&
     );
 
 
@@ -134,7 +134,7 @@ namespace fvm
     tmp<fvMatrix<Type> > laplacian
     (
         const tmp<GeometricField<GType, fvPatchField, volMesh> >&,
-        GeometricField<Type, fvPatchField, volMesh>&,
+        const GeometricField<Type, fvPatchField, volMesh>&,
         const word&
     );
 
@@ -142,7 +142,7 @@ namespace fvm
     tmp<fvMatrix<Type> > laplacian
     (
         const tmp<GeometricField<GType, fvPatchField, volMesh> >&,
-        GeometricField<Type, fvPatchField, volMesh>&
+        const GeometricField<Type, fvPatchField, volMesh>&
     );
 
 
@@ -150,7 +150,7 @@ namespace fvm
     tmp<fvMatrix<Type> > laplacian
     (
         const GeometricField<GType, fvsPatchField, surfaceMesh>&,
-        GeometricField<Type, fvPatchField, volMesh>&,
+        const GeometricField<Type, fvPatchField, volMesh>&,
         const word&
     );
 
@@ -158,7 +158,7 @@ namespace fvm
     tmp<fvMatrix<Type> > laplacian
     (
         const tmp<GeometricField<GType, fvsPatchField, surfaceMesh> >&,
-        GeometricField<Type, fvPatchField, volMesh>&,
+        const GeometricField<Type, fvPatchField, volMesh>&,
         const word&
     );
 
@@ -166,14 +166,14 @@ namespace fvm
     tmp<fvMatrix<Type> > laplacian
     (
         const GeometricField<GType, fvsPatchField, surfaceMesh>&,
-        GeometricField<Type, fvPatchField, volMesh>&
+        const GeometricField<Type, fvPatchField, volMesh>&
     );
 
     template<class Type, class GType>
     tmp<fvMatrix<Type> > laplacian
     (
         const tmp<GeometricField<GType, fvsPatchField, surfaceMesh> >&,
-        GeometricField<Type, fvPatchField, volMesh>&
+        const GeometricField<Type, fvPatchField, volMesh>&
     );
 }
 
diff --git a/src/finiteVolume/finiteVolume/fvm/fvmSup.C b/src/finiteVolume/finiteVolume/fvm/fvmSup.C
index 67c18d3c59..99e2d14746 100644
--- a/src/finiteVolume/finiteVolume/fvm/fvmSup.C
+++ b/src/finiteVolume/finiteVolume/fvm/fvmSup.C
@@ -35,7 +35,7 @@ Foam::tmp<Foam::fvMatrix<Type> >
 Foam::fvm::Su
 (
     const DimensionedField<Type, volMesh>& su,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     const fvMesh& mesh = vf.mesh();
@@ -60,7 +60,7 @@ Foam::tmp<Foam::fvMatrix<Type> >
 Foam::fvm::Su
 (
     const tmp<DimensionedField<Type, volMesh> >& tsu,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > tfvm = fvm::Su(tsu(), vf);
@@ -73,7 +73,7 @@ Foam::tmp<Foam::fvMatrix<Type> >
 Foam::fvm::Su
 (
     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > tfvm = fvm::Su(tsu(), vf);
@@ -86,7 +86,7 @@ Foam::zeroField
 Foam::fvm::Su
 (
     const zeroField&,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     return zeroField();
@@ -98,7 +98,7 @@ Foam::tmp<Foam::fvMatrix<Type> >
 Foam::fvm::Sp
 (
     const DimensionedField<scalar, volMesh>& sp,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     const fvMesh& mesh = vf.mesh();
@@ -123,7 +123,7 @@ Foam::tmp<Foam::fvMatrix<Type> >
 Foam::fvm::Sp
 (
     const tmp<DimensionedField<scalar, volMesh> >& tsp,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > tfvm = fvm::Sp(tsp(), vf);
@@ -136,7 +136,7 @@ Foam::tmp<Foam::fvMatrix<Type> >
 Foam::fvm::Sp
 (
     const tmp<volScalarField>& tsp,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > tfvm = fvm::Sp(tsp(), vf);
@@ -150,7 +150,7 @@ Foam::tmp<Foam::fvMatrix<Type> >
 Foam::fvm::Sp
 (
     const dimensionedScalar& sp,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     const fvMesh& mesh = vf.mesh();
@@ -175,7 +175,7 @@ Foam::zeroField
 Foam::fvm::Sp
 (
     const zeroField&,
-    GeometricField<Type, fvPatchField, volMesh>&
+    const GeometricField<Type, fvPatchField, volMesh>&
 )
 {
     return zeroField();
@@ -187,7 +187,7 @@ Foam::tmp<Foam::fvMatrix<Type> >
 Foam::fvm::SuSp
 (
     const DimensionedField<scalar, volMesh>& susp,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     const fvMesh& mesh = vf.mesh();
@@ -215,7 +215,7 @@ Foam::tmp<Foam::fvMatrix<Type> >
 Foam::fvm::SuSp
 (
     const tmp<DimensionedField<scalar, volMesh> >& tsusp,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > tfvm = fvm::SuSp(tsusp(), vf);
@@ -228,7 +228,7 @@ Foam::tmp<Foam::fvMatrix<Type> >
 Foam::fvm::SuSp
 (
     const tmp<volScalarField>& tsusp,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<fvMatrix<Type> > tfvm = fvm::SuSp(tsusp(), vf);
@@ -241,7 +241,7 @@ Foam::zeroField
 Foam::fvm::SuSp
 (
     const zeroField&,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     return zeroField();
diff --git a/src/finiteVolume/finiteVolume/fvm/fvmSup.H b/src/finiteVolume/finiteVolume/fvm/fvmSup.H
index 425a57c6e3..a5662c1733 100644
--- a/src/finiteVolume/finiteVolume/fvm/fvmSup.H
+++ b/src/finiteVolume/finiteVolume/fvm/fvmSup.H
@@ -56,28 +56,28 @@ namespace fvm
         tmp<fvMatrix<Type> > Su
         (
             const DimensionedField<Type, volMesh>&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         template<class Type>
         tmp<fvMatrix<Type> > Su
         (
             const tmp<DimensionedField<Type, volMesh> >&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         template<class Type>
         tmp<fvMatrix<Type> > Su
         (
             const tmp<GeometricField<Type, fvPatchField, volMesh> >&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         template<class Type>
         zeroField Su
         (
             const zeroField&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
 
@@ -87,21 +87,21 @@ namespace fvm
         tmp<fvMatrix<Type> > Sp
         (
             const DimensionedField<scalar, volMesh>&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         template<class Type>
         tmp<fvMatrix<Type> > Sp
         (
             const tmp<DimensionedField<scalar, volMesh> >&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         template<class Type>
         tmp<fvMatrix<Type> > Sp
         (
             const tmp<volScalarField>&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
 
@@ -109,7 +109,7 @@ namespace fvm
         tmp<fvMatrix<Type> > Sp
         (
             const dimensionedScalar&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
 
@@ -117,7 +117,7 @@ namespace fvm
         zeroField Sp
         (
             const zeroField&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
 
@@ -127,28 +127,28 @@ namespace fvm
         tmp<fvMatrix<Type> > SuSp
         (
             const DimensionedField<scalar, volMesh>&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         template<class Type>
         tmp<fvMatrix<Type> > SuSp
         (
             const tmp<DimensionedField<scalar, volMesh> >&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         template<class Type>
         tmp<fvMatrix<Type> > SuSp
         (
             const tmp<volScalarField>&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         template<class Type>
         zeroField SuSp
         (
             const zeroField&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 }
 
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/extendedLeastSquaresGrad/extendedLeastSquaresGrad.C b/src/finiteVolume/finiteVolume/gradSchemes/extendedLeastSquaresGrad/extendedLeastSquaresGrad.C
index 24b7e8bae0..797e0f5b77 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/extendedLeastSquaresGrad/extendedLeastSquaresGrad.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/extendedLeastSquaresGrad/extendedLeastSquaresGrad.C
@@ -21,7 +21,7 @@ License
     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 "extendedLeastSquaresGrad.H"
@@ -35,27 +35,20 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace fv
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 template<class Type>
-tmp
+Foam::tmp
 <
-    GeometricField
+    Foam::GeometricField
     <
-        typename outerProduct<vector, Type>::type, fvPatchField, volMesh
+        typename Foam::outerProduct<Foam::vector, Type>::type,
+        Foam::fvPatchField,
+    Foam::volMesh
     >
 >
-extendedLeastSquaresGrad<Type>::grad
+Foam::fv::extendedLeastSquaresGrad<Type>::calcGrad
 (
-    const GeometricField<Type, fvPatchField, volMesh>& vsf
+    const GeometricField<Type, fvPatchField, volMesh>& vsf,
+    const word& name
 ) const
 {
     typedef typename outerProduct<vector, Type>::type GradType;
@@ -68,7 +61,7 @@ extendedLeastSquaresGrad<Type>::grad
         (
             IOobject
             (
-                "grad("+vsf.name()+')',
+                name,
                 vsf.instance(),
                 mesh,
                 IOobject::NO_READ,
@@ -120,7 +113,7 @@ extendedLeastSquaresGrad<Type>::grad
 
         if (vsf.boundaryField()[patchi].coupled())
         {
-            Field<Type> neiVsf = 
+            Field<Type> neiVsf =
                 vsf.boundaryField()[patchi].patchNeighbourField();
 
             forAll(neiVsf, patchFaceI)
@@ -162,12 +155,4 @@ extendedLeastSquaresGrad<Type>::grad
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace fv
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/extendedLeastSquaresGrad/extendedLeastSquaresGrad.H b/src/finiteVolume/finiteVolume/gradSchemes/extendedLeastSquaresGrad/extendedLeastSquaresGrad.H
index e8da1ca92f..40003aa703 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/extendedLeastSquaresGrad/extendedLeastSquaresGrad.H
+++ b/src/finiteVolume/finiteVolume/gradSchemes/extendedLeastSquaresGrad/extendedLeastSquaresGrad.H
@@ -102,13 +102,16 @@ public:
 
     // Member Functions
 
-        tmp
+        //- Return the gradient of the given field to the gradScheme::grad
+        //  for optional caching
+        virtual tmp
         <
             GeometricField
             <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
-        > grad
+        > calcGrad
         (
-            const GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>& vsf,
+            const word& name
         ) const;
 };
 
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/fourthGrad/fourthGrad.C b/src/finiteVolume/finiteVolume/gradSchemes/fourthGrad/fourthGrad.C
index f8a4992bd1..c72525e881 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/fourthGrad/fourthGrad.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/fourthGrad/fourthGrad.C
@@ -21,7 +21,7 @@ License
     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 "fourthGrad.H"
@@ -35,27 +35,20 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace fv
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 template<class Type>
-tmp
+Foam::tmp
 <
-    GeometricField
+    Foam::GeometricField
     <
-        typename outerProduct<vector, Type>::type, fvPatchField, volMesh
+        typename Foam::outerProduct<Foam::vector, Type>::type,
+        Foam::fvPatchField,
+        Foam::volMesh
     >
 >
-fourthGrad<Type>::grad
+Foam::fv::fourthGrad<Type>::calcGrad
 (
-    const GeometricField<Type, fvPatchField, volMesh>& vsf
+    const GeometricField<Type, fvPatchField, volMesh>& vsf,
+    const word& name
 ) const
 {
     // The fourth-order gradient is calculated in two passes.  First,
@@ -80,7 +73,7 @@ fourthGrad<Type>::grad
         (
             IOobject
             (
-                "grad("+vsf.name()+')',
+                name,
                 vsf.instance(),
                 mesh,
                 IOobject::NO_READ,
@@ -130,7 +123,7 @@ fourthGrad<Type>::grad
             const scalarField& lambdap = lambda.boundaryField()[patchi];
 
             // Build the d-vectors
-            vectorField pd = 
+            vectorField pd =
                 mesh.Sf().boundaryField()[patchi]
                /(
                    mesh.magSf().boundaryField()[patchi]
@@ -171,12 +164,4 @@ fourthGrad<Type>::grad
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace fv
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/fourthGrad/fourthGrad.H b/src/finiteVolume/finiteVolume/gradSchemes/fourthGrad/fourthGrad.H
index 72e2d3b5ad..1698269b3e 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/fourthGrad/fourthGrad.H
+++ b/src/finiteVolume/finiteVolume/gradSchemes/fourthGrad/fourthGrad.H
@@ -83,13 +83,16 @@ public:
 
     // Member Functions
 
-        tmp
+        //- Return the gradient of the given field to the gradScheme::grad
+        //  for optional caching
+        virtual tmp
         <
             GeometricField
             <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
-        > grad
+        > calcGrad
         (
-            const GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>& vsf,
+            const word& name
         ) const;
 };
 
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C b/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C
index 57dc262370..b8e06f4c55 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C
@@ -21,7 +21,7 @@ License
     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 "gaussGrad.H"
@@ -29,27 +29,20 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace fv
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 template<class Type>
-tmp
+Foam::tmp
 <
-    GeometricField
+    Foam::GeometricField
     <
-        typename outerProduct<vector, Type>::type, fvPatchField, volMesh
+        typename Foam::outerProduct<Foam::vector, Type>::type,
+        Foam::fvPatchField,
+        Foam::volMesh
     >
 >
-gaussGrad<Type>::grad
+Foam::fv::gaussGrad<Type>::gradf
 (
-    const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf,
+    const word& name
 )
 {
     typedef typename outerProduct<vector, Type>::type GradType;
@@ -62,7 +55,7 @@ gaussGrad<Type>::grad
         (
             IOobject
             (
-                "grad("+ssf.name()+')',
+                name,
                 ssf.instance(),
                 mesh,
                 IOobject::NO_READ,
@@ -119,27 +112,29 @@ gaussGrad<Type>::grad
 
 
 template<class Type>
-tmp
+Foam::tmp
 <
-    GeometricField
+    Foam::GeometricField
     <
-        typename outerProduct<vector, Type>::type, fvPatchField, volMesh
+        typename Foam::outerProduct<Foam::vector, Type>::type,
+        Foam::fvPatchField,
+        Foam::volMesh
     >
 >
-gaussGrad<Type>::grad
+Foam::fv::gaussGrad<Type>::calcGrad
 (
-    const GeometricField<Type, fvPatchField, volMesh>& vsf
+    const GeometricField<Type, fvPatchField, volMesh>& vsf,
+    const word& name
 ) const
 {
     typedef typename outerProduct<vector, Type>::type GradType;
 
     tmp<GeometricField<GradType, fvPatchField, volMesh> > tgGrad
     (
-        grad(tinterpScheme_().interpolate(vsf))
+        gradf(tinterpScheme_().interpolate(vsf), name)
     );
     GeometricField<GradType, fvPatchField, volMesh>& gGrad = tgGrad();
 
-    gGrad.rename("grad(" + vsf.name() + ')');
     correctBoundaryConditions(vsf, gGrad);
 
     return tgGrad;
@@ -147,7 +142,7 @@ gaussGrad<Type>::grad
 
 
 template<class Type>
-void gaussGrad<Type>::correctBoundaryConditions
+void Foam::fv::gaussGrad<Type>::correctBoundaryConditions
 (
     const GeometricField<Type, fvPatchField, volMesh>& vsf,
     GeometricField
@@ -174,12 +169,4 @@ void gaussGrad<Type>::correctBoundaryConditions
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace fv
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.H b/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.H
index 8ddfe1bf7c..4b03d582f0 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.H
+++ b/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.H
@@ -89,7 +89,7 @@ public:
             tinterpScheme_(new linear<Type>(mesh))
         {}
 
-        //- Construct from Istream
+        //- Construct from mesh and Istream
         gaussGrad(const fvMesh& mesh, Istream& is)
         :
             gradScheme<Type>(mesh),
@@ -116,31 +116,31 @@ public:
 
     // Member Functions
 
-        //- Return the gradient of the given field 
+        //- Return the gradient of the given field
         //  calculated using Gauss' theorem on the given surface field
         static
         tmp
         <
             GeometricField
             <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
-        > grad
+        > gradf
         (
-            const GeometricField<Type, fvsPatchField, surfaceMesh>&
+            const GeometricField<Type, fvsPatchField, surfaceMesh>&,
+            const word& name
         );
 
-
-        //- Return the gradient of the given field calculated
-        //  using Gauss' theorem on the interpolated field
-        tmp
+        //- Return the gradient of the given field to the gradScheme::grad
+        //  for optional caching
+        virtual tmp
         <
             GeometricField
             <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
-        > grad
+        > calcGrad
         (
-            const GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>& vsf,
+            const word& name
         ) const;
 
-
         //- Correct the boundary values of the gradient using the patchField
         // snGrad functions
         static void correctBoundaryConditions
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrads.C b/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrads.C
index 3aa827351e..2e98e7df2c 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrads.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrads.C
@@ -22,8 +22,6 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-    
 \*---------------------------------------------------------------------------*/
 
 #include "fvMesh.H"
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.C b/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.C
index 6d7cc56560..267896603f 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.C
@@ -22,28 +22,16 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-    Abstract base class for finite volume calculus gradient schemes.
-
 \*---------------------------------------------------------------------------*/
 
 #include "fv.H"
-#include "HashTable.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace fv
-{
+#include "objectRegistry.H"
+#include "solution.H"
 
 // * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
 
 template<class Type>
-tmp<gradScheme<Type> > gradScheme<Type>::New
+Foam::tmp<Foam::fv::gradScheme<Type> > Foam::fv::gradScheme<Type>::New
 (
     const fvMesh& mesh,
     Istream& schemeData
@@ -51,7 +39,8 @@ tmp<gradScheme<Type> > gradScheme<Type>::New
 {
     if (fv::debug)
     {
-        Info<< "gradScheme<Type>::New(Istream& schemeData) : "
+        Info<< "gradScheme<Type>::New"
+               "(const fvMesh& mesh, Istream& schemeData) : "
                "constructing gradScheme<Type>"
             << endl;
     }
@@ -60,7 +49,8 @@ tmp<gradScheme<Type> > gradScheme<Type>::New
     {
         FatalIOErrorIn
         (
-            "gradScheme<Type>::New(Istream& schemeData)",
+            "gradScheme<Type>::New"
+            "(const fvMesh& mesh, Istream& schemeData)",
             schemeData
         )   << "Grad scheme not specified" << endl << endl
             << "Valid grad schemes are :" << endl
@@ -77,7 +67,8 @@ tmp<gradScheme<Type> > gradScheme<Type>::New
     {
         FatalIOErrorIn
         (
-            "gradScheme<Type>::New(Istream& schemeData)",
+            "gradScheme<Type>::New"
+            "(const fvMesh& mesh, Istream& schemeData)",
             schemeData
         )   << "unknown grad scheme " << schemeName << endl << endl
             << "Valid grad schemes are :" << endl
@@ -92,16 +83,151 @@ tmp<gradScheme<Type> > gradScheme<Type>::New
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 template<class Type>
-gradScheme<Type>::~gradScheme()
+Foam::fv::gradScheme<Type>::~gradScheme()
 {}
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace fv
+namespace Foam
+{
+    template<class Type>
+    inline void cachePrintMessage
+    (
+        const char* message,
+        const word& name,
+        const GeometricField<Type, fvPatchField, volMesh>& vf
+    )
+    {
+        if (solution::debug)
+        {
+            Info<< "Cache: " << message << token::SPACE << name
+                << ", " << vf.name() << " event No. " << vf.eventNo()
+                << endl;
+        }
+    }
+}
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+template<class Type>
+Foam::tmp
+<
+    Foam::GeometricField
+    <
+        typename Foam::outerProduct<Foam::vector, Type>::type,
+        Foam::fvPatchField,
+        Foam::volMesh
+    >
+>
+Foam::fv::gradScheme<Type>::grad
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vsf,
+    const word& name
+) const
+{
+    typedef typename outerProduct<vector, Type>::type GradType;
+    typedef GeometricField<GradType, fvPatchField, volMesh> GradFieldType;
+
+    if (!this->mesh().changing() && this->mesh().cache(name))
+    {
+        if (!mesh().objectRegistry::foundObject<GradFieldType>(name))
+        {
+            cachePrintMessage("Calculating and caching", name, vsf);
+            tmp<GradFieldType> tgGrad = calcGrad(vsf, name);
+            regIOobject::store(tgGrad.ptr());
+        }
+
+        cachePrintMessage("Retreiving", name, vsf);
+        GradFieldType& gGrad = const_cast<GradFieldType&>
+        (
+            mesh().objectRegistry::lookupObject<GradFieldType>(name)
+        );
+
+        if (gGrad.upToDate(vsf))
+        {
+            return gGrad;
+        }
+        else
+        {
+            cachePrintMessage("Deleting", name, vsf);
+            gGrad.release();
+            delete &gGrad;
+
+            cachePrintMessage("Recalculating", name, vsf);
+            tmp<GradFieldType> tgGrad = calcGrad(vsf, name);
+
+            cachePrintMessage("Storing", name, vsf);
+            regIOobject::store(tgGrad.ptr());
+            GradFieldType& gGrad = const_cast<GradFieldType&>
+            (
+                mesh().objectRegistry::lookupObject<GradFieldType>(name)
+            );
+
+            return gGrad;
+        }
+    }
+    else
+    {
+        if (mesh().objectRegistry::foundObject<GradFieldType>(name))
+        {
+            GradFieldType& gGrad = const_cast<GradFieldType&>
+            (
+                mesh().objectRegistry::lookupObject<GradFieldType>(name)
+            );
+
+            if (gGrad.ownedByRegistry())
+            {
+                cachePrintMessage("Deleting", name, vsf);
+                gGrad.release();
+                delete &gGrad;
+            }
+        }
+
+        cachePrintMessage("Calculating", name, vsf);
+        return calcGrad(vsf, name);
+    }
+}
+
+
+template<class Type>
+Foam::tmp
+<
+    Foam::GeometricField
+    <
+        typename Foam::outerProduct<Foam::vector, Type>::type,
+        Foam::fvPatchField,
+        Foam::volMesh
+    >
+>
+Foam::fv::gradScheme<Type>::grad
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vsf
+) const
+{
+    return grad(vsf, "grad(" + vsf.name() + ')');
+}
+
+
+template<class Type>
+Foam::tmp
+<
+    Foam::GeometricField
+    <
+        typename Foam::outerProduct<Foam::vector, Type>::type,
+        Foam::fvPatchField,
+        Foam::volMesh
+    >
+>
+Foam::fv::gradScheme<Type>::grad
+(
+    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tvsf
+) const
+{
+    typedef typename outerProduct<vector, Type>::type GradType;
+    typedef GeometricField<GradType, fvPatchField, volMesh> GradFieldType;
+
+    tmp<GradFieldType> tgrad = grad(tvsf());
+    tvsf.clear();
+    return tgrad;
+}
 
-} // End namespace Foam
 
 // ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.H b/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.H
index e338c3dd9e..f1d3fa2262 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.H
+++ b/src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.H
@@ -114,9 +114,8 @@ public:
         );
 
 
-    // Destructor
-
-        virtual ~gradScheme();
+    //- Destructor
+    virtual ~gradScheme();
 
 
     // Member Functions
@@ -127,15 +126,54 @@ public:
             return mesh_;
         }
 
-        //- Calculate and return the grad of the given field
+        //- Calculate and return the grad of the given field.
+        //  Used by grad either to recalculate the cached gradient when it is
+        //  out of date with respect to the field or when it is not cached.
         virtual tmp
+        <
+            GeometricField
+            <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
+        > calcGrad
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&,
+            const word& name
+        ) const = 0;
+
+        //- Calculate and return the grad of the given field
+        //  which may have been cached
+        tmp
+        <
+            GeometricField
+            <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
+        > grad
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&,
+            const word& name
+        ) const;
+
+        //- Calculate and return the grad of the given field
+        //  with the default name
+        //  which may have been cached
+        tmp
         <
             GeometricField
             <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
         > grad
         (
             const GeometricField<Type, fvPatchField, volMesh>&
-        ) const = 0;
+        ) const;
+
+        //- Calculate and return the grad of the given field
+        //  with the default name
+        //  which may have been cached
+        tmp
+        <
+            GeometricField
+            <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
+        > grad
+        (
+            const tmp<GeometricField<Type, fvPatchField, volMesh> >&
+        ) const;
 };
 
 
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.C b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.C
index 2ccdcebfdd..581962315b 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.C
@@ -21,7 +21,7 @@ License
     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 "leastSquaresGrad.H"
@@ -35,27 +35,20 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace fv
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 template<class Type>
-tmp
+Foam::tmp
 <
-    GeometricField
+    Foam::GeometricField
     <
-        typename outerProduct<vector, Type>::type, fvPatchField, volMesh
+        typename Foam::outerProduct<Foam::vector, Type>::type,
+        Foam::fvPatchField,
+        Foam::volMesh
     >
 >
-leastSquaresGrad<Type>::grad
+Foam::fv::leastSquaresGrad<Type>::calcGrad
 (
-    const GeometricField<Type, fvPatchField, volMesh>& vsf
+    const GeometricField<Type, fvPatchField, volMesh>& vsf,
+    const word& name
 ) const
 {
     typedef typename outerProduct<vector, Type>::type GradType;
@@ -68,7 +61,7 @@ leastSquaresGrad<Type>::grad
         (
             IOobject
             (
-                "grad("+vsf.name()+')',
+                name,
                 vsf.instance(),
                 mesh,
                 IOobject::NO_READ,
@@ -116,7 +109,7 @@ leastSquaresGrad<Type>::grad
 
         if (vsf.boundaryField()[patchi].coupled())
         {
-            Field<Type> neiVsf = 
+            Field<Type> neiVsf =
                 vsf.boundaryField()[patchi].patchNeighbourField();
 
             forAll(neiVsf, patchFaceI)
@@ -147,12 +140,4 @@ leastSquaresGrad<Type>::grad
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace fv
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.H b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.H
index 90606b83be..efb15a2465 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.H
+++ b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.H
@@ -89,13 +89,16 @@ public:
 
     // Member Functions
 
-        tmp
+        //- Return the gradient of the given field to the gradScheme::grad
+        //  for optional caching
+        virtual tmp
         <
             GeometricField
             <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
-        > grad
+        > calcGrad
         (
-            const GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>& vsf,
+            const word& name
         ) const;
 };
 
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrad.H b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrad.H
index 4599902885..22083de571 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrad.H
+++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrad.H
@@ -117,18 +117,82 @@ public:
             const Type& extrapolate
         );
 
-
-        tmp
+        //- Return the gradient of the given field to the gradScheme::grad
+        //  for optional caching
+        virtual tmp
         <
             GeometricField
             <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
-        > grad
+        > calcGrad
         (
-            const GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>& vsf,
+            const word& name
         ) const;
 };
 
 
+// * * * * * * * * * * * * Inline Member Function  * * * * * * * * * * * * * //
+
+template<>
+inline void cellLimitedGrad<scalar>::limitFace
+(
+    scalar& limiter,
+    const scalar& maxDelta,
+    const scalar& minDelta,
+    const scalar& extrapolate
+)
+{
+    if (extrapolate > maxDelta + VSMALL)
+    {
+        limiter = min(limiter, maxDelta/extrapolate);
+    }
+    else if (extrapolate < minDelta - VSMALL)
+    {
+        limiter = min(limiter, minDelta/extrapolate);
+    }
+}
+
+
+template<class Type>
+inline void cellLimitedGrad<Type>::limitFace
+(
+    Type& limiter,
+    const Type& maxDelta,
+    const Type& minDelta,
+    const Type& extrapolate
+)
+{
+    for(direction cmpt=0; cmpt<Type::nComponents; cmpt++)
+    {
+        cellLimitedGrad<scalar>::limitFace
+        (
+            limiter.component(cmpt),
+            maxDelta.component(cmpt),
+            minDelta.component(cmpt),
+            extrapolate.component(cmpt)
+        );
+    }
+}
+
+
+// * * * * * * * * Template Member Function Specialisations  * * * * * * * * //
+
+template<>
+tmp<volVectorField> cellLimitedGrad<scalar>::calcGrad
+(
+    const volScalarField& vsf,
+    const word& name
+) const;
+
+
+template<>
+tmp<volTensorField> cellLimitedGrad<vector>::calcGrad
+(
+    const volVectorField& vsf,
+    const word& name
+) const;
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace fv
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrads.C b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrads.C
index c4f0bf524a..fbc2b8fba5 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrads.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrads.C
@@ -36,70 +36,25 @@ License
 
 namespace Foam
 {
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 namespace fv
 {
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-makeFvGradScheme(cellLimitedGrad)
+    makeFvGradScheme(cellLimitedGrad)
+}
+}
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 template<>
-inline void cellLimitedGrad<scalar>::limitFace
+Foam::tmp<Foam::volVectorField>
+Foam::fv::cellLimitedGrad<Foam::scalar>::calcGrad
 (
-    scalar& limiter,
-    const scalar& maxDelta,
-    const scalar& minDelta,
-    const scalar& extrapolate
-)
-{
-    if (extrapolate > maxDelta + VSMALL)
-    {
-        limiter = min(limiter, maxDelta/extrapolate);
-    }
-    else if (extrapolate < minDelta - VSMALL)
-    {
-        limiter = min(limiter, minDelta/extrapolate);
-    }
-}
-
-template<class Type>
-inline void cellLimitedGrad<Type>::limitFace
-(
-    Type& limiter,
-    const Type& maxDelta,
-    const Type& minDelta,
-    const Type& extrapolate
-)
-{
-    for(direction cmpt=0; cmpt<Type::nComponents; cmpt++)
-    {
-        cellLimitedGrad<scalar>::limitFace
-        (
-            limiter.component(cmpt),
-            maxDelta.component(cmpt),
-            minDelta.component(cmpt),
-            extrapolate.component(cmpt)
-        );
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-template<>
-tmp<volVectorField> cellLimitedGrad<scalar>::grad
-(
-    const volScalarField& vsf
+    const volScalarField& vsf,
+    const word& name
 ) const
 {
     const fvMesh& mesh = vsf.mesh();
 
-    tmp<volVectorField> tGrad = basicGradScheme_().grad(vsf);
+    tmp<volVectorField> tGrad = basicGradScheme_().calcGrad(vsf, name);
 
     if (k_ < SMALL)
     {
@@ -244,14 +199,16 @@ tmp<volVectorField> cellLimitedGrad<scalar>::grad
 
 
 template<>
-tmp<volTensorField> cellLimitedGrad<vector>::grad
+Foam::tmp<Foam::volTensorField>
+Foam::fv::cellLimitedGrad<Foam::vector>::calcGrad
 (
-    const volVectorField& vsf
+    const volVectorField& vsf,
+    const word& name
 ) const
 {
     const fvMesh& mesh = vsf.mesh();
 
-    tmp<volTensorField> tGrad = basicGradScheme_().grad(vsf);
+    tmp<volTensorField> tGrad = basicGradScheme_().calcGrad(vsf, name);
 
     if (k_ < SMALL)
     {
@@ -402,12 +359,4 @@ tmp<volTensorField> cellLimitedGrad<vector>::grad
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace fv
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellMDLimitedGrad/cellMDLimitedGrad.H b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellMDLimitedGrad/cellMDLimitedGrad.H
index 05953638b2..d9b389d0f5 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellMDLimitedGrad/cellMDLimitedGrad.H
+++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellMDLimitedGrad/cellMDLimitedGrad.H
@@ -117,13 +117,16 @@ public:
             const vector& dcf
         );
 
-        tmp
+        //- Return the gradient of the given field to the gradScheme::grad
+        //  for optional caching
+        virtual tmp
         <
             GeometricField
             <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
-        > grad
+        > calcGrad
         (
-            const GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>& vsf,
+            const word& name
         ) const;
 };
 
@@ -178,6 +181,24 @@ inline void cellMDLimitedGrad<Type>::limitFace
 }
 
 
+// * * * * * * * * Template Member Function Specialisations  * * * * * * * * //
+
+template<>
+tmp<volVectorField> cellMDLimitedGrad<scalar>::calcGrad
+(
+    const volScalarField& vsf,
+    const word& name
+) const;
+
+
+template<>
+tmp<volTensorField> cellMDLimitedGrad<vector>::calcGrad
+(
+    const volVectorField& vsf,
+    const word& name
+) const;
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace fv
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellMDLimitedGrad/cellMDLimitedGrads.C b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellMDLimitedGrad/cellMDLimitedGrads.C
index 5a35033428..fffc74e97a 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellMDLimitedGrad/cellMDLimitedGrads.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellMDLimitedGrad/cellMDLimitedGrads.C
@@ -36,27 +36,26 @@ License
 
 namespace Foam
 {
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 namespace fv
 {
+    makeFvGradScheme(cellMDLimitedGrad)
+}
+}
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-makeFvGradScheme(cellMDLimitedGrad)
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 template<>
-tmp<volVectorField> cellMDLimitedGrad<scalar>::grad
+Foam::tmp<Foam::volVectorField>
+Foam::fv::cellMDLimitedGrad<Foam::scalar>::calcGrad
 (
-    const volScalarField& vsf
+    const volScalarField& vsf,
+    const word& name
 ) const
 {
     const fvMesh& mesh = vsf.mesh();
 
-    tmp<volVectorField> tGrad = basicGradScheme_().grad(vsf);
+    tmp<volVectorField> tGrad = basicGradScheme_().calcGrad(vsf, name);
 
     if (k_ < SMALL)
     {
@@ -190,14 +189,16 @@ tmp<volVectorField> cellMDLimitedGrad<scalar>::grad
 
 
 template<>
-tmp<volTensorField> cellMDLimitedGrad<vector>::grad
+Foam::tmp<Foam::volTensorField>
+Foam::fv::cellMDLimitedGrad<Foam::vector>::calcGrad
 (
-    const volVectorField& vsf
+    const volVectorField& vsf,
+    const word& name
 ) const
 {
     const fvMesh& mesh = vsf.mesh();
 
-    tmp<volTensorField> tGrad = basicGradScheme_().grad(vsf);
+    tmp<volTensorField> tGrad = basicGradScheme_().calcGrad(vsf, name);
 
     if (k_ < SMALL)
     {
@@ -329,12 +330,4 @@ tmp<volTensorField> cellMDLimitedGrad<vector>::grad
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace fv
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceLimitedGrad/faceLimitedGrad.H b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceLimitedGrad/faceLimitedGrad.H
index 8af5298f42..96393c771e 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceLimitedGrad/faceLimitedGrad.H
+++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceLimitedGrad/faceLimitedGrad.H
@@ -118,17 +118,63 @@ public:
 
     // Member Functions
 
-        tmp
+        //- Return the gradient of the given field to the gradScheme::grad
+        //  for optional caching
+        virtual tmp
         <
             GeometricField
             <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
-        > grad
+        > calcGrad
         (
-            const GeometricField<Type, fvPatchField, volMesh>&
-        ) const;
+            const GeometricField<Type, fvPatchField, volMesh>& vsf,
+            const word& name
+        ) const
+        {
+            return grad(vsf);
+        }
 };
 
 
+// * * * * * * * * * * * * Inline Member Function  * * * * * * * * * * * * * //
+
+template<class Type>
+inline void faceLimitedGrad<Type>::limitFace
+(
+    scalar& limiter,
+    const scalar maxDelta,
+    const scalar minDelta,
+    const scalar extrapolate
+) const
+{
+    if (extrapolate > maxDelta + VSMALL)
+    {
+        limiter = min(limiter, maxDelta/extrapolate);
+    }
+    else if (extrapolate < minDelta - VSMALL)
+    {
+        limiter = min(limiter, minDelta/extrapolate);
+    }
+}
+
+
+// * * * * * * * * Template Member Function Specialisations  * * * * * * * * //
+
+template<>
+tmp<volVectorField> faceLimitedGrad<scalar>::calcGrad
+(
+    const volScalarField& vsf,
+    const word& name
+) const;
+
+
+template<>
+tmp<volTensorField> faceLimitedGrad<vector>::calcGrad
+(
+    const volVectorField& vsf,
+    const word& name
+) const;
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace fv
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceLimitedGrad/faceLimitedGrads.C b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceLimitedGrad/faceLimitedGrads.C
index afaf398fbc..83c065fa54 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceLimitedGrad/faceLimitedGrads.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceLimitedGrad/faceLimitedGrads.C
@@ -21,7 +21,7 @@ License
     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 "faceLimitedGrad.H"
@@ -36,49 +36,26 @@ License
 
 namespace Foam
 {
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 namespace fv
 {
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-makeFvGradScheme(faceLimitedGrad)
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-template<class Type>
-inline void faceLimitedGrad<Type>::limitFace
-(
-    scalar& limiter,
-    const scalar maxDelta,
-    const scalar minDelta,
-    const scalar extrapolate
-) const
-{
-    if (extrapolate > maxDelta + VSMALL)
-    {
-        limiter = min(limiter, maxDelta/extrapolate);
-    }
-    else if (extrapolate < minDelta - VSMALL)
-    {
-        limiter = min(limiter, minDelta/extrapolate);
-    }
+    makeFvGradScheme(faceLimitedGrad)
+}
 }
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 template<>
-tmp<volVectorField> faceLimitedGrad<scalar>::grad
+Foam::tmp<Foam::volVectorField>
+Foam::fv::faceLimitedGrad<Foam::scalar>::calcGrad
 (
-    const volScalarField& vsf
+    const volScalarField& vsf,
+    const word& name
 ) const
 {
     const fvMesh& mesh = vsf.mesh();
 
-    tmp<volVectorField> tGrad = basicGradScheme_().grad(vsf);
+    tmp<volVectorField> tGrad = basicGradScheme_().calcGrad(vsf, name);
 
     if (k_ < SMALL)
     {
@@ -205,14 +182,16 @@ tmp<volVectorField> faceLimitedGrad<scalar>::grad
 
 
 template<>
-tmp<volTensorField> faceLimitedGrad<vector>::grad
+Foam::tmp<Foam::volTensorField>
+Foam::fv::faceLimitedGrad<Foam::vector>::calcGrad
 (
-    const volVectorField& vvf
+    const volVectorField& vvf,
+    const word& name
 ) const
 {
     const fvMesh& mesh = vvf.mesh();
 
-    tmp<volTensorField> tGrad = basicGradScheme_().grad(vvf);
+    tmp<volTensorField> tGrad = basicGradScheme_().calcGrad(vvf, name);
 
     if (k_ < SMALL)
     {
@@ -363,12 +342,4 @@ tmp<volTensorField> faceLimitedGrad<vector>::grad
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace fv
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceMDLimitedGrad/faceMDLimitedGrad.H b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceMDLimitedGrad/faceMDLimitedGrad.H
index f9c826452a..20177019d7 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceMDLimitedGrad/faceMDLimitedGrad.H
+++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceMDLimitedGrad/faceMDLimitedGrad.H
@@ -118,17 +118,38 @@ public:
 
     // Member Functions
 
-        tmp
+        //- Return the gradient of the given field to the gradScheme::grad
+        //  for optional caching
+        virtual tmp
         <
             GeometricField
             <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
-        > grad
+        > calcGrad
         (
-            const GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>& vsf,
+            const word& name
         ) const;
 };
 
 
+// * * * * * * * * Template Member Function Specialisations  * * * * * * * * //
+
+template<>
+tmp<volVectorField> faceMDLimitedGrad<scalar>::calcGrad
+(
+    const volScalarField& vsf,
+    const word& name
+) const;
+
+
+template<>
+tmp<volTensorField> faceMDLimitedGrad<vector>::calcGrad
+(
+    const volVectorField& vsf,
+    const word& name
+) const;
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace fv
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceMDLimitedGrad/faceMDLimitedGrads.C b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceMDLimitedGrad/faceMDLimitedGrads.C
index 27124c4673..87ddd25096 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceMDLimitedGrad/faceMDLimitedGrads.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceMDLimitedGrad/faceMDLimitedGrads.C
@@ -21,7 +21,7 @@ License
     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 "faceMDLimitedGrad.H"
@@ -37,28 +37,25 @@ License
 
 namespace Foam
 {
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 namespace fv
 {
+    makeFvGradScheme(faceMDLimitedGrad)
+}
+}
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-makeFvGradScheme(faceMDLimitedGrad)
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-// FaceLimited scalar gradient
 template<>
-tmp<volVectorField> faceMDLimitedGrad<scalar>::grad
+Foam::tmp<Foam::volVectorField>
+Foam::fv::faceMDLimitedGrad<Foam::scalar>::calcGrad
 (
-    const volScalarField& vsf
+    const volScalarField& vsf,
+    const word& name
 ) const
 {
     const fvMesh& mesh = vsf.mesh();
 
-    tmp<volVectorField> tGrad = basicGradScheme_().grad(vsf);
+    tmp<volVectorField> tGrad = basicGradScheme_().calcGrad(vsf, name);
 
     if (k_ < SMALL)
     {
@@ -189,14 +186,16 @@ tmp<volVectorField> faceMDLimitedGrad<scalar>::grad
 
 
 template<>
-tmp<volTensorField> faceMDLimitedGrad<vector>::grad
+Foam::tmp<Foam::volTensorField>
+Foam::fv::faceMDLimitedGrad<Foam::vector>::calcGrad
 (
-    const volVectorField& vvf
+    const volVectorField& vvf,
+    const word& name
 ) const
 {
     const fvMesh& mesh = vvf.mesh();
 
-    tmp<volTensorField> tGrad = basicGradScheme_().grad(vvf);
+    tmp<volTensorField> tGrad = basicGradScheme_().calcGrad(vvf, name);
 
     if (k_ < SMALL)
     {
@@ -327,12 +326,4 @@ tmp<volTensorField> faceMDLimitedGrad<vector>::grad
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace fv
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianScheme.C b/src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianScheme.C
index 1df55b7e17..b58c2bd41d 100644
--- a/src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianScheme.C
+++ b/src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianScheme.C
@@ -47,7 +47,7 @@ tmp<fvMatrix<Type> >
 gaussLaplacianScheme<Type, GType>::fvmLaplacianUncorrected
 (
     const surfaceScalarField& gammaMagSf,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     tmp<surfaceScalarField> tdeltaCoeffs =
@@ -67,9 +67,9 @@ gaussLaplacianScheme<Type, GType>::fvmLaplacianUncorrected
     fvm.upper() = deltaCoeffs.internalField()*gammaMagSf.internalField();
     fvm.negSumDiag();
 
-    forAll(fvm.psi().boundaryField(), patchI)
+    forAll(vf.boundaryField(), patchI)
     {
-        const fvPatchField<Type>& psf = fvm.psi().boundaryField()[patchI];
+        const fvPatchField<Type>& psf = vf.boundaryField()[patchI];
         const fvsPatchScalarField& patchGamma =
             gammaMagSf.boundaryField()[patchI];
 
@@ -149,7 +149,7 @@ tmp<fvMatrix<Type> >
 gaussLaplacianScheme<Type, GType>::fvmLaplacian
 (
     const GeometricField<GType, fvsPatchField, surfaceMesh>& gamma,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     const fvMesh& mesh = this->mesh();
diff --git a/src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianScheme.H b/src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianScheme.H
index a6c3138833..12b33ae0e4 100644
--- a/src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianScheme.H
+++ b/src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianScheme.H
@@ -62,7 +62,7 @@ class gaussLaplacianScheme
         tmp<fvMatrix<Type> > fvmLaplacianUncorrected
         (
             const surfaceScalarField& gammaMagSf,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > gammaSnGradCorr
@@ -126,7 +126,7 @@ public:
         tmp<fvMatrix<Type> > fvmLaplacian
         (
             const GeometricField<GType, fvsPatchField, surfaceMesh>&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         tmp<GeometricField<Type, fvPatchField, volMesh> > fvcLaplacian
@@ -146,7 +146,7 @@ template<>                                                                  \
 tmp<fvMatrix<Type> > gaussLaplacianScheme<Type, scalar>::fvmLaplacian       \
 (                                                                           \
     const GeometricField<scalar, fvsPatchField, surfaceMesh>&,              \
-    GeometricField<Type, fvPatchField, volMesh>&                            \
+    const GeometricField<Type, fvPatchField, volMesh>&                      \
 );                                                                          \
                                                                             \
 template<>                                                                  \
diff --git a/src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianSchemes.C b/src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianSchemes.C
index 5dab5ceef0..782f14ea4a 100644
--- a/src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianSchemes.C
+++ b/src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianSchemes.C
@@ -21,7 +21,7 @@ License
     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 "gaussLaplacianScheme.H"
@@ -44,7 +44,7 @@ Foam::tmp<Foam::fvMatrix<Foam::Type> >                                       \
 Foam::fv::gaussLaplacianScheme<Foam::Type, Foam::scalar>::fvmLaplacian       \
 (                                                                            \
     const GeometricField<scalar, fvsPatchField, surfaceMesh>& gamma,         \
-    GeometricField<Type, fvPatchField, volMesh>& vf                          \
+    const GeometricField<Type, fvPatchField, volMesh>& vf                    \
 )                                                                            \
 {                                                                            \
     const fvMesh& mesh = this->mesh();                                       \
diff --git a/src/finiteVolume/finiteVolume/laplacianSchemes/laplacianScheme/laplacianScheme.C b/src/finiteVolume/finiteVolume/laplacianSchemes/laplacianScheme/laplacianScheme.C
index eb14f253ad..78ce0bc2c6 100644
--- a/src/finiteVolume/finiteVolume/laplacianSchemes/laplacianScheme/laplacianScheme.C
+++ b/src/finiteVolume/finiteVolume/laplacianSchemes/laplacianScheme/laplacianScheme.C
@@ -102,7 +102,7 @@ tmp<fvMatrix<Type> >
 laplacianScheme<Type, GType>::fvmLaplacian
 (
     const GeometricField<GType, fvPatchField, volMesh>& gamma,
-    GeometricField<Type, fvPatchField, volMesh>& vf
+    const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
     return fvmLaplacian(tinterpGammaScheme_().interpolate(gamma)(), vf);
diff --git a/src/finiteVolume/finiteVolume/laplacianSchemes/laplacianScheme/laplacianScheme.H b/src/finiteVolume/finiteVolume/laplacianSchemes/laplacianScheme/laplacianScheme.H
index b36d87df9e..4bbe44ed9a 100644
--- a/src/finiteVolume/finiteVolume/laplacianSchemes/laplacianScheme/laplacianScheme.H
+++ b/src/finiteVolume/finiteVolume/laplacianSchemes/laplacianScheme/laplacianScheme.H
@@ -173,13 +173,13 @@ public:
         virtual tmp<fvMatrix<Type> > fvmLaplacian
         (
             const GeometricField<GType, fvsPatchField, surfaceMesh>&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         ) = 0;
 
         virtual tmp<fvMatrix<Type> > fvmLaplacian
         (
             const GeometricField<GType, fvPatchField, volMesh>&,
-            GeometricField<Type, fvPatchField, volMesh>&
+            const GeometricField<Type, fvPatchField, volMesh>&
         );
 
         virtual tmp<GeometricField<Type, fvPatchField, volMesh> > fvcLaplacian
diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrad.C b/src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrad.C
index 0168bd626a..631d6a9c2f 100644
--- a/src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrad.C
+++ b/src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrad.C
@@ -22,9 +22,6 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-    Simple central-difference snGrad scheme with non-orthogonal correction.
-
 \*---------------------------------------------------------------------------*/
 
 #include "correctedSnGrad.H"
@@ -34,28 +31,44 @@ Description
 #include "fvcGrad.H"
 #include "gaussGrad.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace fv
-{
-
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 template<class Type>
-correctedSnGrad<Type>::~correctedSnGrad()
+Foam::fv::correctedSnGrad<Type>::~correctedSnGrad()
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
-tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
-correctedSnGrad<Type>::correction
+Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> >
+Foam::fv::correctedSnGrad<Type>::fullGradCorrection
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+) const
+{
+    const fvMesh& mesh = this->mesh();
+
+    // construct GeometricField<Type, fvsPatchField, surfaceMesh>
+    tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tssf =
+        mesh.correctionVectors()
+      & linear<typename outerProduct<vector, Type>::type>(mesh).interpolate
+        (
+            gradScheme<Type>::New
+            (
+                mesh,
+                mesh.gradScheme(vf.name())
+            )().grad(vf, "grad(" + vf.name() + ')')
+        );
+    tssf().rename("snGradCorr(" + vf.name() + ')');
+
+    return tssf;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> >
+Foam::fv::correctedSnGrad<Type>::correction
 (
     const GeometricField<Type, fvPatchField, volMesh>& vf
 ) const
@@ -89,7 +102,7 @@ correctedSnGrad<Type>::correction
             mesh.correctionVectors()
           & linear
             <
-                typename 
+                typename
                 outerProduct<vector, typename pTraits<Type>::cmptType>::type
             >(mesh).interpolate
             (
@@ -108,12 +121,4 @@ correctedSnGrad<Type>::correction
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace fv
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrad.H b/src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrad.H
index 2358867066..3eaa42bfdf 100644
--- a/src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrad.H
+++ b/src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrad.H
@@ -108,12 +108,36 @@ public:
         }
 
         //- Return the explicit correction to the correctedSnGrad
-        //  for the given field
+        //  for the given field using the gradient of the field
+        tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
+        fullGradCorrection
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const;
+
+        //- Return the explicit correction to the correctedSnGrad
+        //  for the given field using the gradients of the field components
         virtual tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
         correction(const GeometricField<Type, fvPatchField, volMesh>&) const;
 };
 
 
+// * * * * * * * * Template Member Function Specialisations  * * * * * * * * //
+
+template<>
+tmp<surfaceScalarField> correctedSnGrad<scalar>::correction
+(
+    const volScalarField& vsf
+) const;
+
+
+template<>
+tmp<surfaceVectorField> correctedSnGrad<vector>::correction
+(
+    const volVectorField& vvf
+) const;
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace fv
diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrads.C b/src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrads.C
index b416a806fa..8377ba7ef4 100644
--- a/src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrads.C
+++ b/src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrads.C
@@ -22,9 +22,6 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-    Simple central-difference snGrad scheme with non-orthogonal correction.
-
 \*---------------------------------------------------------------------------*/
 
 #include "correctedSnGrad.H"
@@ -40,4 +37,27 @@ namespace fv
 }
 }
 
+
+template<>
+Foam::tmp<Foam::surfaceScalarField>
+Foam::fv::correctedSnGrad<Foam::scalar>::correction
+(
+    const volScalarField& vsf
+) const
+{
+    return fullGradCorrection(vsf);
+}
+
+
+template<>
+Foam::tmp<Foam::surfaceVectorField>
+Foam::fv::correctedSnGrad<Foam::vector>::correction
+(
+    const volVectorField& vvf
+) const
+{
+    return fullGradCorrection(vvf);
+}
+
+
 // ************************************************************************* //
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
index 1e49de4683..74b4282638 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
@@ -183,7 +183,7 @@ void Foam::fvMatrix<Type>::addBoundarySource
 template<class Type>
 Foam::fvMatrix<Type>::fvMatrix
 (
-    GeometricField<Type, fvPatchField, volMesh>& psi,
+    const GeometricField<Type, fvPatchField, volMesh>& psi,
     const dimensionSet& ds
 )
 :
@@ -227,7 +227,13 @@ Foam::fvMatrix<Type>::fvMatrix
         );
     }
 
-    psi_.boundaryField().updateCoeffs();
+    // Update the boundary coefficients of psi without changing it's event No.
+    GeometricField<Type, fvPatchField, volMesh>& psiRef =
+       const_cast<GeometricField<Type, fvPatchField, volMesh>&>(psi_);
+
+    label currentStatePsi = psiRef.eventNo();
+    psiRef.boundaryField().updateCoeffs();
+    psiRef.eventNo() = currentStatePsi;
 }
 
 
@@ -322,7 +328,7 @@ Foam::fvMatrix<Type>::fvMatrix(const tmp<fvMatrix<Type> >& tfvm)
 template<class Type>
 Foam::fvMatrix<Type>::fvMatrix
 (
-    GeometricField<Type, fvPatchField, volMesh>& psi,
+    const GeometricField<Type, fvPatchField, volMesh>& psi,
     Istream& is
 )
 :
@@ -404,12 +410,17 @@ void Foam::fvMatrix<Type>::setValues
     const unallocLabelList& nei = mesh.neighbour();
 
     scalarField& Diag = diag();
+    Field<Type>& psi =
+        const_cast
+        <
+            GeometricField<Type, fvPatchField, volMesh>&
+        >(psi_).internalField();
 
     forAll(cellLabels, i)
     {
         label celli = cellLabels[i];
 
-        psi_[celli] = values[i];
+        psi[celli] = values[i];
         source_[celli] = values[i]*Diag[celli];
 
         if (symmetric() || asymmetric())
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
index 4fbc2430d7..8569862e13 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
@@ -117,8 +117,9 @@ public:
 
     // Private data
 
-        // Reference to GeometricField<Type, fvPatchField, volMesh>
-        GeometricField<Type, fvPatchField, volMesh>& psi_;
+        //- Const reference to GeometricField<Type, fvPatchField, volMesh>
+        //  Converted into a non-const reference at the point of solution.
+        const GeometricField<Type, fvPatchField, volMesh>& psi_;
 
         //- Dimension set
         dimensionSet dimensions_;
@@ -237,7 +238,7 @@ public:
         //- Construct given a field to solve for
         fvMatrix
         (
-            GeometricField<Type, fvPatchField, volMesh>&,
+            const GeometricField<Type, fvPatchField, volMesh>&,
             const dimensionSet&
         );
 
@@ -245,12 +246,12 @@ public:
         fvMatrix(const fvMatrix<Type>&);
 
         //- Construct as copy of tmp<fvMatrix<Type> > deleting argument
-#       ifdef ConstructFromTmp
+        #ifdef ConstructFromTmp
         fvMatrix(const tmp<fvMatrix<Type> >&);
-#       endif
+        #endif
 
         //- Construct from Istream given field to solve for
-        fvMatrix(GeometricField<Type, fvPatchField, volMesh>&, Istream&);
+        fvMatrix(const GeometricField<Type, fvPatchField, volMesh>&, Istream&);
 
 
     // Destructor
@@ -267,11 +268,6 @@ public:
                 return psi_;
             }
 
-            GeometricField<Type, fvPatchField, volMesh>& psi()
-            {
-                return psi_;
-            }
-
             const dimensionSet& dimensions() const
             {
                 return dimensions_;
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C
index 21f3738a02..c56899bd62 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C
@@ -63,10 +63,13 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve
             << endl;
     }
 
+    GeometricField<Type, fvPatchField, volMesh>& psi =
+       const_cast<GeometricField<Type, fvPatchField, volMesh>&>(psi_);
+
     lduMatrix::solverPerformance solverPerfVec
     (
         "fvMatrix<Type>::solve",
-        psi_.name()
+        psi.name()
     );
 
     scalarField saveDiag = diag();
@@ -82,7 +85,7 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve
     (
         pow
         (
-            psi_.mesh().solutionD(),
+            psi.mesh().solutionD(),
             pTraits<typename powProduct<Vector<label>, Type::rank>::type>::zero
         )
     );
@@ -93,7 +96,7 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve
 
         // copy field and source
 
-        scalarField psiCmpt = psi_.internalField().component(cmpt);
+        scalarField psiCmpt = psi.internalField().component(cmpt);
         addBoundaryDiag(diag(), cmpt);
 
         scalarField sourceCmpt = source.component(cmpt);
@@ -109,7 +112,7 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve
         );
 
         lduInterfaceFieldPtrsList interfaces =
-            psi_.boundaryField().interfaces();
+            psi.boundaryField().interfaces();
 
         // Use the initMatrixInterfaces and updateMatrixInterfaces to correct
         // bouCoeffsCmpt for the explicit part of the coupled boundary
@@ -137,7 +140,7 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve
         // Solver call
         solverPerf = lduMatrix::solver::New
         (
-            psi_.name() + pTraits<Type>::componentNames[cmpt],
+            psi.name() + pTraits<Type>::componentNames[cmpt],
             *this,
             bouCoeffsCmpt,
             intCoeffsCmpt,
@@ -156,11 +159,11 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve
             solverPerfVec = solverPerf;
         }
 
-        psi_.internalField().replace(cmpt, psiCmpt);
+        psi.internalField().replace(cmpt, psiCmpt);
         diag() = saveDiag;
     }
 
-    psi_.correctBoundaryConditions();
+    psi.correctBoundaryConditions();
 
     return solverPerfVec;
 }
diff --git a/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C b/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C
index 35e5c3ac35..5ea0ba6a17 100644
--- a/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C
+++ b/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C
@@ -99,6 +99,10 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Foam::scalar>::fvSolver::solve
     const dictionary& solverControls
 )
 {
+    GeometricField<scalar, fvPatchField, volMesh>& psi =
+        const_cast<GeometricField<scalar, fvPatchField, volMesh>&>
+        (fvMat_.psi());
+
     scalarField saveDiag = fvMat_.diag();
     fvMat_.addBoundaryDiag(fvMat_.diag(), 0);
 
@@ -108,14 +112,17 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Foam::scalar>::fvSolver::solve
     // assign new solver controls
     solver_->read(solverControls);
 
-    lduMatrix::solverPerformance solverPerf =
-        solver_->solve(fvMat_.psi().internalField(), totalSource);
+    lduMatrix::solverPerformance solverPerf = solver_->solve
+    (
+        psi.internalField(),
+        totalSource
+    );
 
     solverPerf.print();
 
     fvMat_.diag() = saveDiag;
 
-    fvMat_.psi().correctBoundaryConditions();
+    psi.correctBoundaryConditions();
 
     return solverPerf;
 }
@@ -134,6 +141,9 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Foam::scalar>::solve
             << endl;
     }
 
+    GeometricField<scalar, fvPatchField, volMesh>& psi =
+       const_cast<GeometricField<scalar, fvPatchField, volMesh>&>(psi_);
+
     scalarField saveDiag = diag();
     addBoundaryDiag(diag(), 0);
 
@@ -143,19 +153,19 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Foam::scalar>::solve
     // Solver call
     lduMatrix::solverPerformance solverPerf = lduMatrix::solver::New
     (
-        psi_.name(),
+        psi.name(),
         *this,
         boundaryCoeffs_,
         internalCoeffs_,
-        psi_.boundaryField().interfaces(),
+        psi.boundaryField().interfaces(),
         solverControls
-    )->solve(psi_.internalField(), totalSource);
+    )->solve(psi.internalField(), totalSource);
 
     solverPerf.print();
 
     diag() = saveDiag;
 
-    psi_.correctBoundaryConditions();
+    psi.correctBoundaryConditions();
 
     return solverPerf;
 }
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/LimitedScheme.C b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/LimitedScheme.C
index 1465c2648e..ff2ecdbdc3 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/LimitedScheme.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/LimitedScheme.C
@@ -66,8 +66,10 @@ tmp<surfaceScalarField> LimitedScheme<Type, Limiter, LimitFunc>::limiter
     const GeometricField<typename Limiter::phiType, fvPatchField, volMesh>&
         lPhi = tlPhi();
 
-    GeometricField<typename Limiter::gradPhiType, fvPatchField, volMesh>
-        gradc(fvc::grad(lPhi));
+    tmp<GeometricField<typename Limiter::gradPhiType, fvPatchField, volMesh> >
+        tgradc(fvc::grad(lPhi));
+    const GeometricField<typename Limiter::gradPhiType, fvPatchField, volMesh>&
+        gradc = tgradc();
 
     const surfaceScalarField& CDweights = mesh.surfaceInterpolation::weights();
 
@@ -116,7 +118,7 @@ tmp<surfaceScalarField> LimitedScheme<Type, Limiter, LimitFunc>::limiter
                 gradc.boundaryField()[patchi].patchNeighbourField();
 
             // Build the d-vectors
-            vectorField pd = 
+            vectorField pd =
                 mesh.Sf().boundaryField()[patchi]
                /(
                    mesh.magSf().boundaryField()[patchi]
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/linearUpwind/linearUpwind.C b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/linearUpwind/linearUpwind.C
index 4cda19451b..104cbdd29b 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/linearUpwind/linearUpwind.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/linearUpwind/linearUpwind.C
@@ -67,9 +67,22 @@ Foam::linearUpwind<Type>::correction
     const volVectorField& C = mesh.C();
     const surfaceVectorField& Cf = mesh.Cf();
 
-    GeometricField
-        <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
-        gradVf = gradScheme_().grad(vf);
+    tmp
+    <
+        GeometricField
+        <
+            typename outerProduct<vector, Type>::type,
+            fvPatchField,
+            volMesh
+        >
+    > tgradVf = gradScheme_().grad(vf, gradSchemeName_);
+
+    const GeometricField
+    <
+        typename outerProduct<vector, Type>::type,
+        fvPatchField,
+        volMesh
+    >& gradVf = tgradVf();
 
     forAll(faceFlux, facei)
     {
@@ -95,7 +108,7 @@ Foam::linearUpwind<Type>::correction
 
         if (pSfCorr.coupled())
         {
-            const unallocLabelList& pOwner = 
+            const unallocLabelList& pOwner =
                 mesh.boundary()[patchi].faceCells();
 
             const vectorField& pCf = Cf.boundaryField()[patchi];
@@ -106,7 +119,7 @@ Foam::linearUpwind<Type>::correction
                 gradVf.boundaryField()[patchi].patchNeighbourField();
 
             // Build the d-vectors
-            vectorField pd = 
+            vectorField pd =
                 mesh.Sf().boundaryField()[patchi]
                /(
                    mesh.magSf().boundaryField()[patchi]
@@ -129,7 +142,7 @@ Foam::linearUpwind<Type>::correction
                 }
                 else
                 {
-                    pSfCorr[facei] = 
+                    pSfCorr[facei] =
                         (pCf[facei] - pd[facei] - C[own]) & pGradVfNei[facei];
                 }
             }
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/linearUpwind/linearUpwind.H b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/linearUpwind/linearUpwind.H
index 2b51590f85..e5e278a7e9 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/linearUpwind/linearUpwind.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/linearUpwind/linearUpwind.H
@@ -56,6 +56,7 @@ class linearUpwind
 {
     // Private Data
 
+        word gradSchemeName_;
         tmp<fv::gradScheme<Type> > gradScheme_;
 
 
@@ -84,6 +85,7 @@ public:
         )
         :
             upwind<Type>(mesh, faceFlux),
+            gradSchemeName_("grad"),
             gradScheme_
             (
                 new fv::gaussGrad<Type>(mesh)
@@ -100,12 +102,13 @@ public:
         )
         :
             upwind<Type>(mesh, schemeData),
+            gradSchemeName_(schemeData),
             gradScheme_
             (
                 fv::gradScheme<Type>::New
                 (
                     mesh,
-                    schemeData
+                    mesh.gradScheme(gradSchemeName_)
                 )
             )
         {}
@@ -119,12 +122,13 @@ public:
         )
         :
             upwind<Type>(mesh, faceFlux, schemeData),
+            gradSchemeName_(schemeData),
             gradScheme_
             (
                 fv::gradScheme<Type>::New
                 (
                     mesh,
-                    schemeData
+                    mesh.gradScheme(gradSchemeName_)
                 )
             )
         {}
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/linearUpwind/linearUpwindV.C b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/linearUpwind/linearUpwindV.C
index 80a137b396..c4b52f4da0 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/linearUpwind/linearUpwindV.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/linearUpwind/linearUpwindV.C
@@ -74,9 +74,22 @@ Foam::linearUpwindV<Type>::correction
     const vectorField& C = mesh.C();
     const vectorField& Cf = mesh.Cf();
 
-    GeometricField
-        <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
-        gradVf = gradScheme_().grad(vf);
+    tmp
+    <
+        GeometricField
+        <
+            typename outerProduct<vector, Type>::type,
+            fvPatchField,
+            volMesh
+        >
+    > tgradVf = gradScheme_().grad(vf, gradSchemeName_);
+
+    const GeometricField
+    <
+        typename outerProduct<vector, Type>::type,
+        fvPatchField,
+        volMesh
+    >& gradVf = tgradVf();
 
     forAll(faceFlux, facei)
     {
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/linearUpwind/linearUpwindV.H b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/linearUpwind/linearUpwindV.H
index 7f8c9de742..9d477610a7 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/linearUpwind/linearUpwindV.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/linearUpwind/linearUpwindV.H
@@ -56,6 +56,7 @@ class linearUpwindV
 {
     // Private Data
 
+        word gradSchemeName_;
         tmp<fv::gradScheme<Type> > gradScheme_;
 
 
@@ -84,6 +85,7 @@ public:
         )
         :
             upwind<Type>(mesh, faceFlux),
+            gradSchemeName_("grad"),
             gradScheme_
             (
                 new fv::gaussGrad<Type>(mesh)
@@ -100,12 +102,13 @@ public:
         )
         :
             upwind<Type>(mesh, schemeData),
+            gradSchemeName_(schemeData),
             gradScheme_
             (
                 fv::gradScheme<Type>::New
                 (
                     mesh,
-                    schemeData
+                    mesh.gradScheme(gradSchemeName_)
                 )
             )
         {}
@@ -119,12 +122,13 @@ public:
         )
         :
             upwind<Type>(mesh, faceFlux, schemeData),
+            gradSchemeName_(schemeData),
             gradScheme_
             (
                 fv::gradScheme<Type>::New
                 (
                     mesh,
-                    schemeData
+                    mesh.gradScheme(gradSchemeName_)
                 )
             )
         {}
diff --git a/src/lagrangian/intermediate/Make/files b/src/lagrangian/intermediate/Make/files
index 270bc40e2e..6695181469 100644
--- a/src/lagrangian/intermediate/Make/files
+++ b/src/lagrangian/intermediate/Make/files
@@ -47,6 +47,7 @@ $(REACTINGMPPARCEL)/makeBasicReactingMultiphaseParcelSubmodels.C
 /* bolt-on models */
 submodels/addOns/radiation/absorptionEmission/cloudAbsorptionEmission/cloudAbsorptionEmission.C
 submodels/addOns/radiation/scatter/cloudScatter/cloudScatter.C
+submodels/Kinematic/PatchInteractionModel/LocalInteraction/patchInteractionData.C
 
 
 /* data entries */
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index 1ae536e4db..8c020d3739 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -286,7 +286,13 @@ bool Foam::KinematicParcel<ParcelType>::hitPatch
     ParcelType& p = static_cast<ParcelType&>(*this);
     td.cloud().postProcessing().postPatch(p, patchI);
 
-    return td.cloud().patchInteraction().correct(pp, this->face(), U_);
+    return td.cloud().patchInteraction().correct
+    (
+        pp,
+        this->face(),
+        td.keepParticle,
+        U_
+    );
 }
 
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/LocalInteraction.C b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/LocalInteraction.C
index ffe1f09631..3beaa8abaf 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/LocalInteraction.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/LocalInteraction.C
@@ -26,7 +26,7 @@ License
 
 #include "LocalInteraction.H"
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
 template <class CloudType>
 Foam::label Foam::LocalInteraction<CloudType>::applyToPatch
@@ -46,7 +46,7 @@ Foam::label Foam::LocalInteraction<CloudType>::applyToPatch
 }
 
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * * //
 
 template <class CloudType>
 Foam::LocalInteraction<CloudType>::LocalInteraction
@@ -62,6 +62,7 @@ Foam::LocalInteraction<CloudType>::LocalInteraction
     const polyMesh& mesh = cloud.mesh();
     const polyBoundaryMesh& bMesh = mesh.boundaryMesh();
 
+    // check that user patches are valid region patches
     forAll(patchData_, patchI)
     {
         const word& patchName = patchData_[patchI].patchName();
@@ -70,7 +71,7 @@ Foam::LocalInteraction<CloudType>::LocalInteraction
         {
             FatalErrorIn("LocalInteraction(const dictionary&, CloudType&)")
                 << "Patch " << patchName << " not found. Available patches "
-                << "are: " << bMesh.names() << exit(FatalError);
+                << "are: " << bMesh.names() << nl << exit(FatalError);
         }
     }
 
@@ -95,6 +96,26 @@ Foam::LocalInteraction<CloudType>::LocalInteraction
             << "interaction. Please specify data for patches:" << nl
             << badWalls << nl << exit(FatalError);
     }
+
+    // check that interactions are valid/specified
+    forAll(patchData_, patchI)
+    {
+        const word& interactionTypeName =
+            patchData_[patchI].interactionTypeName();
+        const typename PatchInteractionModel<CloudType>::interactionType& it =
+            this->wordToInteractionType(interactionTypeName);
+
+        if (it == PatchInteractionModel<CloudType>::itOther)
+        {
+            const word& patchName = patchData_[patchI].patchName();
+            FatalErrorIn("LocalInteraction(const dictionary&, CloudType&)")
+                << "Unknown patch interaction type "
+                << interactionTypeName << " for patch " << patchName
+                << ". Valid selections are:"
+                << this->PatchInteractionModel<CloudType>::interactionTypeNames_
+                << nl << exit(FatalError);
+        }
+    }
 }
 
 
@@ -105,7 +126,7 @@ Foam::LocalInteraction<CloudType>::~LocalInteraction()
 {}
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 template<class CloudType>
 bool Foam::LocalInteraction<CloudType>::active() const
@@ -119,6 +140,7 @@ bool Foam::LocalInteraction<CloudType>::correct
 (
     const polyPatch& pp,
     const label faceId,
+    bool& keepParticle,
     vector& U
 ) const
 {
@@ -126,18 +148,64 @@ bool Foam::LocalInteraction<CloudType>::correct
 
     if (patchI >= 0)
     {
-        vector nw = pp.faceAreas()[pp.whichFace(faceId)];
-        nw /= mag(nw);
+        typename PatchInteractionModel<CloudType>::interactionType it =
+            this->wordToInteractionType
+            (
+                patchData_[patchI].interactionTypeName()
+            );
 
-        scalar Un = U & nw;
-        vector Ut = U - Un*nw;
-
-        if (Un > 0)
+        switch (it)
         {
-            U -= (1.0 + patchData_[patchI].e())*Un*nw;
-        }
+            case PatchInteractionModel<CloudType>::itEscape:
+            {
+                keepParticle = false;
+                U = vector::zero;
+                break;
+            }
+            case PatchInteractionModel<CloudType>::itStick:
+            {
+                keepParticle = true;
+                U = vector::zero;
+                break;
+            }
+            case PatchInteractionModel<CloudType>::itRebound:
+            {
+                keepParticle = true;
 
-        U -= patchData_[patchI].mu()*Ut;
+                vector nw = pp.faceAreas()[pp.whichFace(faceId)];
+                nw /= mag(nw);
+
+                scalar Un = U & nw;
+                vector Ut = U - Un*nw;
+
+                if (Un > 0)
+                {
+                    U -= (1.0 + patchData_[patchI].e())*Un*nw;
+                }
+
+                U -= patchData_[patchI].mu()*Ut;
+
+                break;
+            }
+            default:
+            {
+                FatalErrorIn
+                (
+                    "bool LocalInteraction<CloudType>::correct"
+                    "("
+                        "const polyPatch&, "
+                        "const label, "
+                        "bool&, "
+                        "vector&"
+                    ") const"
+                )   << "Unknown interaction type "
+                    << patchData_[patchI].interactionTypeName()
+                    << "(" << it << ") for patch "
+                    << patchData_[patchI].patchName()
+                    << ". Valid selections are:" << this->interactionTypeNames_
+                    << endl << abort(FatalError);
+            }
+        }
 
         return true;
     }
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/LocalInteraction.H b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/LocalInteraction.H
index 49d36e9c6d..902c4316e5 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/LocalInteraction.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/LocalInteraction.H
@@ -34,7 +34,7 @@ Description
 #define LocalInteraction_H
 
 #include "PatchInteractionModel.H"
-#include "dictionaryEntry.H"
+#include "patchInteractionData.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -49,78 +49,6 @@ class LocalInteraction
 :
     public PatchInteractionModel<CloudType>
 {
-    class patchInteractionData
-    {
-        // Private data
-
-            //- Patch name
-            word patchName_;
-
-            //- Elasticity coefficient
-            scalar e_;
-
-            //- Restitution coefficient
-            scalar mu_;
-
-
-    public:
-
-        //- Construct null
-        patchInteractionData()
-        :
-            patchName_("unknownPatch"),
-            e_(0.0),
-            mu_(0.0)
-        {}
-
-        //- Construct from dictionary
-        patchInteractionData(const dictionary& dict);
-
-        // Member functions
-
-            // Access
-
-                //- Return const access to the patch name
-                const word& patchName() const
-                {
-                    return patchName_;
-                }
-
-                //- Return const access to the elasticity coefficient
-                scalar e() const
-                {
-                    return e_;
-                }
-
-                //- Return const access to the restitution coefficient
-                scalar mu() const
-                {
-                    return mu_;
-                }
-
-
-            // I-O
-
-                //- Istream operator
-                friend Istream& operator>>(Istream& is, patchInteractionData& pid)
-                {
-                    is.check
-                    (
-                        "Istream& operator>>"
-                        "(Istream&, patchInteractionData&)"
-                    );
-
-                    const dictionaryEntry entry(dictionary::null, is);
-
-                    pid.patchName_ = entry.keyword();
-                    entry.lookup("e") >> pid.e_;
-                    entry.lookup("mu") >> pid.mu_;
-
-                    return is;
-                }
-    };
-
-
     // Private data
 
         //- List of participating patches
@@ -164,6 +92,7 @@ public:
         (
             const polyPatch& pp,
             const label faceId,
+            bool& keepParticle,
             vector& U
         ) const;
 };
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/patchInteractionData.C b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/patchInteractionData.C
new file mode 100644
index 0000000000..f5569fa94d
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/patchInteractionData.C
@@ -0,0 +1,89 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-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 "patchInteractionData.H"
+#include "dictionaryEntry.H"
+#include "PatchInteractionModel.H"
+
+// * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * * //
+
+Foam::patchInteractionData::patchInteractionData()
+:
+    interactionTypeName_("unknownInteractionTypeName"),
+    patchName_("unknownPatch"),
+    e_(0.0),
+    mu_(0.0)
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+const Foam::word& Foam::patchInteractionData::interactionTypeName() const
+{
+    return interactionTypeName_;
+}
+
+
+const Foam::word& Foam::patchInteractionData::patchName() const
+{
+    return patchName_;
+}
+
+
+Foam::scalar Foam::patchInteractionData::e() const
+{
+    return e_;
+}
+
+
+Foam::scalar Foam::patchInteractionData::mu() const
+{
+    return mu_;
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Istream& Foam::operator>>
+(
+    Istream& is,
+    patchInteractionData& pid
+)
+{
+    is.check("Istream& operator>>(Istream&, patchInteractionData&)");
+
+    const dictionaryEntry entry(dictionary::null, is);
+
+    pid.patchName_ = entry.keyword();
+    entry.lookup("type") >> pid.interactionTypeName_;
+    pid.e_ = entry.lookupOrDefault<scalar>("e", 1.0);
+    pid.mu_ = entry.lookupOrDefault<scalar>("mu", 0.0);
+
+    return is;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/patchInteractionData.H b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/patchInteractionData.H
new file mode 100644
index 0000000000..d53bb3dcf6
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/patchInteractionData.H
@@ -0,0 +1,157 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-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::patchInteractionData
+
+Description
+    Helper class for the LocalInteraction patch interaction model
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef patchInteractionData_H
+#define patchInteractionData_H
+
+#include "Istream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                    Class patchInteractionData Declaration
+\*---------------------------------------------------------------------------*/
+
+// Forward declaration of classes
+class patchInteractionData;
+
+// Forward declaration of friend functions
+Istream& operator>>
+(
+    Istream& is,
+    patchInteractionData& pid
+);
+
+
+class patchInteractionData
+{
+    // Private data
+
+        //- Interaction type name
+        word interactionTypeName_;
+
+        //- Patch name
+        word patchName_;
+
+        //- Elasticity coefficient
+        scalar e_;
+
+        //- Restitution coefficient
+        scalar mu_;
+
+
+public:
+
+    // Constructor
+
+        //- Construct null
+        patchInteractionData();
+
+
+    // Member functions
+
+        // Access
+
+            //- Return const access to the interaction type name
+            const word& interactionTypeName() const;
+
+            //- Return const access to the patch name
+            const word& patchName() const;
+
+            //- Return const access to the elasticity coefficient
+            scalar e() const;
+
+            //- Return const access to the restitution coefficient
+            scalar mu() const;
+
+
+        // I-O
+
+            //- Istream operator
+            friend Istream& operator>>
+            (
+                Istream& is,
+                patchInteractionData& pid
+            );
+/*            {
+                is.check
+                (
+                    "Istream& operator>>"
+                    "(Istream&, patchInteractionData&)"
+                );
+
+                const dictionaryEntry entry(dictionary::null, is);
+
+                pid.patchName_ = entry.keyword();
+                entry.lookup("type") >> pid.interactionTypeName_;
+                pid.e_ = entry.lookupOrDefault<scalar>("e", 1.0);
+                pid.mu_ = entry.lookupOrDefault<scalar>("mu", 0.0);
+
+                if
+                (
+                    PatchInteractionModel<CloudType>::wordToInteractionType
+                    (
+                        pid.interactionTypeName_
+                    )
+                == PatchInteractionModel<CloudType>::itOther)
+                {
+                    FatalErrorIn
+                    (
+                        "friend Istream& operator>>"
+                        "("
+                            "Istream&, "
+                            "patchInteractionData&"
+                        ")"
+                    )   << "Unknown patch interaction type "
+                        << pid.interactionTypeName_
+                        << ". Valid selections are:"
+                        << PatchInteractionModel<CloudType>::
+                        interactionTypeNames_
+                        << endl << abort(FatalError);
+                }
+
+                return is;
+            }
+*/};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/PatchInteractionModel/PatchInteractionModel.C b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/PatchInteractionModel/PatchInteractionModel.C
index 117e6005c8..873d5fb346 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/PatchInteractionModel/PatchInteractionModel.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/PatchInteractionModel/PatchInteractionModel.C
@@ -26,6 +26,76 @@ License
 
 #include "PatchInteractionModel.H"
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::wordList Foam::PatchInteractionModel<CloudType>::interactionTypeNames_
+(
+    IStringStream
+    (
+        "(rebound stick escape)"
+    )()
+);
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::word Foam::PatchInteractionModel<CloudType>::interactionTypeToWord
+(
+    const interactionType& itEnum
+)
+{
+    switch (itEnum)
+    {
+        case itRebound:
+        {
+            return "rebound";
+            break;
+        }
+        case itStick:
+        {
+            return "stick";
+            break;
+        }
+        case itEscape:
+        {
+            return "escape";
+            break;
+        }
+        default:
+        {
+            return "other";
+        }
+    }
+}
+
+
+template<class CloudType>
+typename Foam::PatchInteractionModel<CloudType>::interactionType
+Foam::PatchInteractionModel<CloudType>::wordToInteractionType
+(
+    const word& itWord
+)
+{
+    if (itWord == "rebound")
+    {
+        return itRebound;
+    }
+    else if (itWord == "stick")
+    {
+        return itStick;
+    }
+    else if (itWord == "escape")
+    {
+        return itEscape;
+    }
+    else
+    {
+        return itOther;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/PatchInteractionModel/PatchInteractionModel.H b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/PatchInteractionModel/PatchInteractionModel.H
index 986d30d3a1..cf9eb2b05c 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/PatchInteractionModel/PatchInteractionModel.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/PatchInteractionModel/PatchInteractionModel.H
@@ -40,6 +40,7 @@ SourceFiles
 #include "IOdictionary.H"
 #include "autoPtr.H"
 #include "runTimeSelectionTables.H"
+#include "polyPatch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -53,6 +54,24 @@ namespace Foam
 template<class CloudType>
 class PatchInteractionModel
 {
+public:
+
+    // Public enumerations
+
+        // Interaction types
+        enum interactionType
+        {
+            itRebound,
+            itStick,
+            itEscape,
+            itOther
+        };
+
+        static wordList interactionTypeNames_;
+
+
+private:
+
     // Private data
 
         //- The cloud dictionary
@@ -121,6 +140,12 @@ public:
 
     // Member Functions
 
+        //- Convert interaction result to word
+        static word interactionTypeToWord(const interactionType& itEnum);
+
+        //- Convert word to interaction result
+        static interactionType wordToInteractionType(const word& itWord);
+
         //- Flag to indicate whether model activates patch interaction model
         virtual bool active() const = 0;
 
@@ -130,6 +155,7 @@ public:
         (
             const polyPatch& pp,
             const label faceId,
+            bool& keepParticle,
             vector& U
         ) const = 0;
 };
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/Rebound/Rebound.C b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/Rebound/Rebound.C
index e5cfdd3540..44d8915497 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/Rebound/Rebound.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/Rebound/Rebound.C
@@ -61,9 +61,12 @@ bool Foam::Rebound<CloudType>::correct
 (
     const polyPatch& pp,
     const label faceId,
+    bool& keepParticle,
     vector& U
 ) const
 {
+    keepParticle = true;
+
     vector nw = pp.faceAreas()[pp.whichFace(faceId)];
     nw /= mag(nw);
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/Rebound/Rebound.H b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/Rebound/Rebound.H
index 7e3126996f..3980ff9c4d 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/Rebound/Rebound.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/Rebound/Rebound.H
@@ -82,6 +82,7 @@ public:
         (
             const polyPatch& pp,
             const label faceId,
+            bool& keepParticle,
             vector& U
         ) const;
 };
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/StandardWallInteraction/StandardWallInteraction.C b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/StandardWallInteraction/StandardWallInteraction.C
index f87c0dd08b..fcf8f04a2f 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/StandardWallInteraction/StandardWallInteraction.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/StandardWallInteraction/StandardWallInteraction.C
@@ -36,9 +36,45 @@ Foam::StandardWallInteraction<CloudType>::StandardWallInteraction
 )
 :
     PatchInteractionModel<CloudType>(dict, cloud, typeName),
-    e_(dimensionedScalar(this->coeffDict().lookup("e")).value()),
-    mu_(dimensionedScalar(this->coeffDict().lookup("mu")).value())
-{}
+    interactionType_
+    (
+        this->wordToInteractionType(this->coeffDict().lookup("type"))
+    ),
+    e_(0.0),
+    mu_(0.0)
+{
+    switch (interactionType_)
+    {
+        case PatchInteractionModel<CloudType>::itOther:
+        {
+            word interactionTypeName(this->coeffDict().lookup("type"));
+
+            FatalErrorIn
+            (
+                "StandardWallInteraction<CloudType>::StandardWallInteraction"
+                "("
+                    "const dictionary&, "
+                    "CloudType& cloud"
+                ")"
+            )   << "Unknown interaction result type "
+                << interactionTypeName
+                << ". Valid selections are:" << this->interactionTypeNames_
+                << endl << exit(FatalError);
+
+            break;
+        }
+        case PatchInteractionModel<CloudType>::itRebound:
+        {
+            e_ = this->coeffDict().lookupOrDefault("e", 1.0);
+            mu_ = this->coeffDict().lookupOrDefault("mu", 0.0);
+            break;
+        }
+        default:
+        {
+            // do nothing
+        }
+    }
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -62,23 +98,62 @@ bool Foam::StandardWallInteraction<CloudType>::correct
 (
     const polyPatch& pp,
     const label faceId,
+    bool& keepParticle,
     vector& U
 ) const
 {
     if (isA<wallPolyPatch>(pp))
     {
-        vector nw = pp.faceAreas()[pp.whichFace(faceId)];
-        nw /= mag(nw);
-
-        scalar Un = U & nw;
-        vector Ut = U - Un*nw;
-
-        if (Un > 0)
+        switch (interactionType_)
         {
-            U -= (1.0 + e_)*Un*nw;
-        }
+            case PatchInteractionModel<CloudType>::itEscape:
+            {
+                keepParticle = false;
+                U = vector::zero;
+                break;
+            }
+            case PatchInteractionModel<CloudType>::itStick:
+            {
+                keepParticle = true;
+                U = vector::zero;
+                break;
+            }
+            case PatchInteractionModel<CloudType>::itRebound:
+            {
+                keepParticle = true;
 
-        U -= mu_*Ut;
+                vector nw = pp.faceAreas()[pp.whichFace(faceId)];
+                nw /= mag(nw);
+
+                scalar Un = U & nw;
+                vector Ut = U - Un*nw;
+
+                if (Un > 0)
+                {
+                    U -= (1.0 + e_)*Un*nw;
+                }
+
+                U -= mu_*Ut;
+
+                break;
+            }
+            default:
+            {
+                FatalErrorIn
+                (
+                    "bool StandardWallInteraction<CloudType>::correct"
+                    "("
+                        "const polyPatch&, "
+                        "const label, "
+                        "bool&, "
+                        "vector&"
+                    ") const"
+                )   << "Unknown interaction type "
+                    << this->interactionTypeToWord(interactionType_)
+                    << "(" << interactionType_ << ")" << endl
+                    << abort(FatalError);
+            }
+        }
 
         return true;
     }
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/StandardWallInteraction/StandardWallInteraction.H b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/StandardWallInteraction/StandardWallInteraction.H
index 3e2175e469..ab1f198bdc 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/StandardWallInteraction/StandardWallInteraction.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/StandardWallInteraction/StandardWallInteraction.H
@@ -26,7 +26,19 @@ Class
     Foam::StandardWallInteraction
 
 Description
-    Wall interaction based on restitution and elasticity coefficients
+    Wall interaction model. Three choices:
+    - rebound - optionally specify elasticity and resitution coefficients
+    - stick   - particles assigined zero velocity
+    - escape  - remove particle from the domain
+
+    Example usage:
+
+        StandardWallInteractionCoeffs
+        {
+            type        rebound; // stick, escape
+            e           1;       // optional - elasticity coeff
+            mu          0;       // optional - restitution coeff
+        }
 
 \*---------------------------------------------------------------------------*/
 
@@ -40,7 +52,7 @@ Description
 namespace Foam
 {
 /*---------------------------------------------------------------------------*\
-                      Class StandardWallInteraction Declaration
+                   Class StandardWallInteraction Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class CloudType>
@@ -48,13 +60,19 @@ class StandardWallInteraction
 :
     public PatchInteractionModel<CloudType>
 {
-    // Private data
+protected:
 
-        //- Elasticity
-        const scalar e_;
+    // Protected data
+
+        //- Interaction type
+        typename PatchInteractionModel<CloudType>::interactionType
+            interactionType_;
+
+        //- Elasticity coefficient
+        scalar e_;
 
         //- Restitution coefficient
-        const scalar mu_;
+        scalar mu_;
 
 
 public:
@@ -84,6 +102,7 @@ public:
         (
             const polyPatch& pp,
             const label faceId,
+            bool& keepParticle,
             vector& U
         ) const;
 };
diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
index d763866a49..3927ad7740 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
+++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C
@@ -161,7 +161,7 @@ void Foam::LiquidEvaporation<CloudType>::calculate
         // vapour diffusivity [m2/s]
         scalar Dab = liquids_->properties()[lid].D(pc, Ts);
 
-        // Saturation pressure for species i [pa]
+        // saturation pressure for species i [pa]
         // - carrier phase pressure assumed equal to the liquid vapour pressure
         //   close to the surface
         // NOTE: if pSat > pc then particle is superheated
diff --git a/src/mesh/blockMesh/curvedEdges/BSpline.C b/src/mesh/blockMesh/curvedEdges/BSpline.C
index 5aaf693c29..3e166f25cb 100644
--- a/src/mesh/blockMesh/curvedEdges/BSpline.C
+++ b/src/mesh/blockMesh/curvedEdges/BSpline.C
@@ -46,7 +46,7 @@ Foam::pointField Foam::BSpline::findKnots
     register const scalar oneSixth = 1.0/6.0;
     register const scalar twoThird = 2.0/3.0;
 
-    simpleMatrix<vector> M(NKnots+2);
+    simpleMatrix<vector> M(NKnots+2, 0, vector::zero);
 
     // set up the matrix
     M[0][0] = -0.5*scalar(NKnots - 1);
diff --git a/src/mesh/blockMesh/curvedEdges/spline.C b/src/mesh/blockMesh/curvedEdges/spline.C
index 67b8c79559..cddf6fa8ed 100644
--- a/src/mesh/blockMesh/curvedEdges/spline.C
+++ b/src/mesh/blockMesh/curvedEdges/spline.C
@@ -28,7 +28,7 @@ License
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
-Foam::scalar Foam::spline::B(const scalar& tau)
+Foam::scalar Foam::spline::B(const scalar tau)
 {
     if (tau <= -2.0 || tau >= 2.0)
     {
@@ -52,7 +52,7 @@ Foam::scalar Foam::spline::B(const scalar& tau)
     }
     else
     {
-        FatalErrorIn("spline::B(const scalar&)")
+        FatalErrorIn("spline::B(const scalar)")
             << "Programming error???, "
             << "tau = " << tau
             << abort(FatalError);
diff --git a/src/mesh/blockMesh/curvedEdges/spline.H b/src/mesh/blockMesh/curvedEdges/spline.H
index 8609a64674..74e73afdc6 100644
--- a/src/mesh/blockMesh/curvedEdges/spline.H
+++ b/src/mesh/blockMesh/curvedEdges/spline.H
@@ -60,7 +60,7 @@ class spline
     // Private Member Functions
 
         //- Blending function for constructing spline
-        static scalar B(const scalar&);
+        static scalar B(const scalar);
 
         //- Disallow default bitwise copy construct
         spline(const spline&);
diff --git a/src/meshTools/octree/treeBoundBox.H b/src/meshTools/octree/treeBoundBox.H
index 69791e04a8..04b420a531 100644
--- a/src/meshTools/octree/treeBoundBox.H
+++ b/src/meshTools/octree/treeBoundBox.H
@@ -158,10 +158,6 @@ public:
         static const FixedList<vector, 6> faceNormals;
 
 
-        //- Face on which neighbour is
-        static direction neighbourFaceBits(const label&);
-
-
     // Constructors
 
         //- Construct null setting points to zero
diff --git a/src/surfMesh/MeshedSurface/MeshedSurface.C b/src/surfMesh/MeshedSurface/MeshedSurface.C
index 8bcab42089..1c1622c0eb 100644
--- a/src/surfMesh/MeshedSurface/MeshedSurface.C
+++ b/src/surfMesh/MeshedSurface/MeshedSurface.C
@@ -507,7 +507,7 @@ void Foam::MeshedSurface<Face>::movePoints(const pointField& newPoints)
 
 
 template<class Face>
-void Foam::MeshedSurface<Face>::scalePoints(const scalar& scaleFactor)
+void Foam::MeshedSurface<Face>::scalePoints(const scalar scaleFactor)
 {
     // avoid bad scaling
     if (scaleFactor > 0 && scaleFactor != 1.0)
diff --git a/src/surfMesh/MeshedSurface/MeshedSurface.H b/src/surfMesh/MeshedSurface/MeshedSurface.H
index 5677233e0b..ea7125fa83 100644
--- a/src/surfMesh/MeshedSurface/MeshedSurface.H
+++ b/src/surfMesh/MeshedSurface/MeshedSurface.H
@@ -338,7 +338,7 @@ public:
         virtual void movePoints(const pointField&);
 
         //- Scale points. A non-positive factor is ignored
-        virtual void scalePoints(const scalar&);
+        virtual void scalePoints(const scalar);
 
         //- Reset primitive data (points, faces and zones)
         //  Note, optimized to avoid overwriting data (with Xfer::null)
diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C b/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C
index 95cfc76388..9813c796bb 100644
--- a/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C
+++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C
@@ -63,11 +63,8 @@ Foam::scalar Foam::EulerImplicit<CompType, ThermoType>::solve
     const scalar dt
 ) const
 {
-    scalar pf, cf, pr, cr;
-    label lRef, rRef;
-
-    label nSpecie = this->model_.nSpecie();
-    simpleMatrix<scalar> RR(nSpecie);
+    const label nSpecie = this->model_.nSpecie();
+    simpleMatrix<scalar> RR(nSpecie, 0, 0);
 
     for (label i=0; i<nSpecie; i++)
     {
@@ -83,6 +80,9 @@ Foam::scalar Foam::EulerImplicit<CompType, ThermoType>::solve
     {
         const Reaction<ThermoType>& R = this->model_.reactions()[i];
 
+        scalar pf, cf, pr, cr;
+        label lRef, rRef;
+
         scalar omegai = this->model_.omega
         (
             R, c, T, p, pf, cf, lRef, pr, cr, rRef
@@ -133,7 +133,7 @@ Foam::scalar Foam::EulerImplicit<CompType, ThermoType>::solve
 
     // estimate the next time step
     scalar tMin = GREAT;
-    label nEqns = this->model_.nEqns();
+    const label nEqns = this->model_.nEqns();
     scalarField c1(nEqns, 0.0);
 
     for (label i=0; i<nSpecie; i++)
diff --git a/src/triSurface/triSurface/triSurface.C b/src/triSurface/triSurface/triSurface.C
index 3b8d2f9d5f..65258557e8 100644
--- a/src/triSurface/triSurface/triSurface.C
+++ b/src/triSurface/triSurface/triSurface.C
@@ -33,8 +33,6 @@ License
 #include "SortableList.H"
 #include "PackedBoolList.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 namespace Foam
 {
@@ -783,7 +781,6 @@ const Foam::labelList& Foam::triSurface::edgeOwner() const
 }
 
 
-//- Move points
 void Foam::triSurface::movePoints(const pointField& newPoints)
 {
     // Remove all geometry dependent data
@@ -797,8 +794,7 @@ void Foam::triSurface::movePoints(const pointField& newPoints)
 }
 
 
-// scale points
-void Foam::triSurface::scalePoints(const scalar& scaleFactor)
+void Foam::triSurface::scalePoints(const scalar scaleFactor)
 {
     // avoid bad scaling
     if (scaleFactor > 0 && scaleFactor != 1.0)
diff --git a/src/triSurface/triSurface/triSurface.H b/src/triSurface/triSurface/triSurface.H
index 9ee5b5ca4d..7ef87db1ab 100644
--- a/src/triSurface/triSurface/triSurface.H
+++ b/src/triSurface/triSurface/triSurface.H
@@ -324,7 +324,7 @@ public:
             virtual void movePoints(const pointField&);
 
             //- Scale points. A non-positive factor is ignored
-            virtual void scalePoints(const scalar&);
+            virtual void scalePoints(const scalar);
 
             //- Check/remove duplicate/degenerate triangles
             void checkTriangles(const bool verbose);
diff --git a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/coalCloud1Properties b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/coalCloud1Properties
index f270546201..ef11516391 100644
--- a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/coalCloud1Properties
+++ b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/coalCloud1Properties
@@ -107,8 +107,9 @@ ManualInjectionCoeffs
 
 StandardWallInteractionCoeffs
 {
-    e               e [ 0 0 0 0 0 ] 1;
-    mu              mu [ 0 0 0 0 0 ] 0;
+    type            rebound;
+    e               1;
+    mu              0;
 }
 
 RanzMarshallCoeffs
diff --git a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/limestoneCloud1Properties b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/limestoneCloud1Properties
index a17be1e36f..f4cf89a700 100644
--- a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/limestoneCloud1Properties
+++ b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/limestoneCloud1Properties
@@ -92,8 +92,9 @@ ManualInjectionCoeffs
 
 StandardWallInteractionCoeffs
 {
-    e               e [ 0 0 0 0 0 ] 1;
-    mu              mu [ 0 0 0 0 0 ] 0;
+    type            rebound;
+    e               1;
+    mu              0;
 }
 
 RanzMarshallCoeffs
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/reactingCloud1Properties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/reactingCloud1Properties
index 5c99dbf554..b1b89e17eb 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/reactingCloud1Properties
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/constant/reactingCloud1Properties
@@ -110,8 +110,7 @@ ReactingLookupTableInjectionCoeffs
 
 StandardWallInteractionCoeffs
 {
-    e               e [ 0 0 0 0 0 ] 1;
-    mu              mu [ 0 0 0 0 0 ] 0;
+    type            rebound;
 }
 
 LocalInteractionCoeffs
@@ -120,13 +119,11 @@ LocalInteractionCoeffs
     (
         walls
         {
-            e       1;
-            mu      0;
+            type    rebound;
         }
         cycLeft
         {
-            e       1;
-            mu      0;
+            type    rebound;
         }
     );
 }
diff --git a/tutorials/lagrangian/reactingParcelFoam/evaporationTest/constant/reactingCloud1Properties b/tutorials/lagrangian/reactingParcelFoam/evaporationTest/constant/reactingCloud1Properties
index 95f4df9ceb..97619b95a7 100644
--- a/tutorials/lagrangian/reactingParcelFoam/evaporationTest/constant/reactingCloud1Properties
+++ b/tutorials/lagrangian/reactingParcelFoam/evaporationTest/constant/reactingCloud1Properties
@@ -99,8 +99,7 @@ ManualInjectionCoeffs
 
 StandardWallInteractionCoeffs
 {
-    e               e [ 0 0 0 0 0 ] 1;
-    mu              mu [ 0 0 0 0 0 ] 0;
+    type            rebound;
 }
 
 RanzMarshallCoeffs
diff --git a/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/kinematicCloud1Properties b/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/kinematicCloud1Properties
index f88d7e1090..469aaa10f6 100644
--- a/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/kinematicCloud1Properties
+++ b/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/kinematicCloud1Properties
@@ -104,8 +104,7 @@ ConeInjectionCoeffs
 
 StandardWallInteractionCoeffs
 {
-    e               e [ 0 0 0 0 0 ] 1;
-    mu              mu [ 0 0 0 0 0 ] 0;
+    type            rebound;
 }
 
 
diff --git a/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/thermoCloud1Properties b/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/thermoCloud1Properties
index 67823f6aab..59402e4ea3 100644
--- a/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/thermoCloud1Properties
+++ b/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/constant/thermoCloud1Properties
@@ -92,8 +92,7 @@ ManualInjectionCoeffs
 
 StandardWallInteractionCoeffs
 {
-    e               e [ 0 0 0 0 0 ] 1;
-    mu              mu [ 0 0 0 0 0 ] 0;
+    type            rebound;
 }
 
 RanzMarshallCoeffs