Merge branch 'master' of ssh://noisy/home/noisy3/OpenFOAM/OpenFOAM-dev

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

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>&);

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_)
     {

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

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)
     {

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

src/OpenFOAM/global/unitConversion/unitConversion.H

@@ -45,13 +45,13 @@ namespace Foam
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 //- Conversion from degrees to radians
-inline scalar degToRad(const scalar& deg)
+inline scalar degToRad(const scalar deg)
 {
     return (deg*pi/180.0);
 }
 
 //- Conversion from radians to degrees
-inline scalar radToDeg(const scalar& rad)
+inline scalar radToDeg(const scalar rad)
 {
     return (rad*180.0/pi);
 }

src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C

@@ -32,7 +32,20 @@ template<class Type>
 Foam::simpleMatrix<Type>::simpleMatrix(const label mSize)
 :
     scalarSquareMatrix(mSize),
-    source_(mSize, pTraits<Type>::zero)
+    source_(mSize)
+{}
+
+
+template<class Type>
+Foam::simpleMatrix<Type>::simpleMatrix
+(
+    const label mSize,
+    const scalar coeffVal,
+    const Type& sourceVal
+)
+:
+    scalarSquareMatrix(mSize, mSize, coeffVal),
+    source_(mSize, sourceVal)
 {}
 
 

src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H

@@ -26,7 +26,7 @@ Class
     Foam::simpleMatrix
 
 Description
-    Foam::simpleMatrix
+    A simple square matrix solver with scalar coefficients.
 
 SourceFiles
     simpleMatrix.C
@@ -75,8 +75,13 @@ public:
     // Constructors
 
         //- Construct given size
+        //  Note: this does not initialise the coefficients or the source.
         simpleMatrix(const label);
 
+        //- Construct given size and initial values for the
+        //  coefficients and source
+        simpleMatrix(const label, const scalar, const Type&);
+
         //- Construct from components
         simpleMatrix(const scalarSquareMatrix&, const Field<Type>&);
 
@@ -91,11 +96,13 @@ public:
 
         // Access
 
+            //- Return access to the source
             Field<Type>& source()
             {
                 return source_;
             }
 
+            //- Return const-access to the source
             const Field<Type>& source() const
             {
                 return source_;

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"))

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

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)

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

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)
     {

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

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)
     {

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

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

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

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"

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

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

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

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

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

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>

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

src/finiteVolume/finiteVolume/d2dt2Schemes/EulerD2dt2Scheme/EulerD2dt2Scheme.C

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

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

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

src/finiteVolume/finiteVolume/d2dt2Schemes/steadyStateD2dt2Scheme/steadyStateD2dt2Scheme.C

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

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

src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C

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

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;

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 =

src/finiteVolume/finiteVolume/ddtSchemes/CrankNicholsonDdtScheme/CrankNicholsonDdtScheme.H

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

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

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;

src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C

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

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;

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

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;

src/finiteVolume/finiteVolume/ddtSchemes/boundedBackwardDdtScheme/boundedBackwardDdtScheme.C

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

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

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;
 
 

src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.C

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

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;

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

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

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

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

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

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

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

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

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

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

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

src/finiteVolume/finiteVolume/gradSchemes/extendedLeastSquaresGrad/extendedLeastSquaresGrad.C

@@ -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,
@@ -162,12 +155,4 @@ extendedLeastSquaresGrad<Type>::grad
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace fv
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //

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

src/finiteVolume/finiteVolume/gradSchemes/fourthGrad/fourthGrad.C

@@ -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,
@@ -171,12 +164,4 @@ fourthGrad<Type>::grad
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace fv
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //

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

src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C

@@ -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
-
 // ************************************************************************* //

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),
@@ -123,24 +123,24 @@ public:
         <
             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 to the gradScheme::grad
-        //- Return the gradient of the given field calculated
+        //  for optional caching
-        //  using Gauss' theorem on the interpolated field
+        virtual tmp
-        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

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"

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"
+#include "objectRegistry.H"
-
+#include "solution.H"
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace fv
-{
 
 // * * * * * * * * * * * * * * * * * 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
 
 // ************************************************************************* //

src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.H

@@ -114,9 +114,8 @@ public:
         );
 
 
-    // Destructor
+    //- Destructor
-
+    virtual ~gradScheme();
-        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;
 };
 
 

src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.C

@@ -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,
@@ -147,12 +140,4 @@ leastSquaresGrad<Type>::grad
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace fv
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //

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

src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrad.H

@@ -117,18 +117,82 @@ public:
             const Type& extrapolate
         );
 
-
+        //- Return the gradient of the given field to the gradScheme::grad
-        tmp
+        //  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

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 volScalarField& vsf,
-    const scalar& maxDelta,
+    const word& name
-    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
 {
     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
-
 // ************************************************************************* //

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

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
-
 // ************************************************************************* //

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 GeometricField<Type, fvPatchField, volMesh>& vsf,
-        ) const;
+            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

src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceLimitedGrad/faceLimitedGrads.C

@@ -36,49 +36,26 @@ License
 
 namespace Foam
 {
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 namespace fv
 {
-
+    makeFvGradScheme(faceLimitedGrad)
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+}
-
-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);
-    }
 }
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 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
-
 // ************************************************************************* //

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

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

@@ -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
-
 // ************************************************************************* //

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

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

src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianSchemes.C

@@ -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();                                       \

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

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

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> >
+Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> >
-correctedSnGrad<Type>::correction
+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
@@ -108,12 +121,4 @@ correctedSnGrad<Type>::correction
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace fv
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //

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

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);
+}
+
+
 // ************************************************************************* //

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

src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H

@@ -117,8 +117,9 @@ public:
 
     // Private data
 
-        // Reference to GeometricField<Type, fvPatchField, volMesh>
+        //- Const reference to GeometricField<Type, fvPatchField, volMesh>
-        GeometricField<Type, fvPatchField, volMesh>& psi_;
+        //  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_;

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

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 =
+    lduMatrix::solverPerformance solverPerf = solver_->solve
-        solver_->solve(fvMat_.psi().internalField(), totalSource);
+    (
+        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;
 }

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>
+    tmp<GeometricField<typename Limiter::gradPhiType, fvPatchField, volMesh> >
-        gradc(fvc::grad(lPhi));
+        tgradc(fvc::grad(lPhi));
+    const GeometricField<typename Limiter::gradPhiType, fvPatchField, volMesh>&
+        gradc = tgradc();
 
     const surfaceScalarField& CDweights = mesh.surfaceInterpolation::weights();
 

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
+    tmp
-        <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
+    <
-        gradVf = gradScheme_().grad(vf);
+        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)
     {

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_)
                 )
             )
         {}

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
+    tmp
-        <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
+    <
-        gradVf = gradScheme_().grad(vf);
+        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)
     {

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_)
                 )
             )
         {}

src/mesh/blockMesh/curvedEdges/BSpline.C

@@ -48,7 +48,7 @@ Foam::pointField Foam::BSpline::findKnots
     register scalar oneSixth = 1.0/6.0;
     register scalar twoThird = 2.0/3.0;
 
-    simpleMatrix<vector> M(newnKnots);
+    simpleMatrix<vector> M(newnKnots, 0, vector::zero);
 
     // set up the matrix
 

src/mesh/blockMesh/curvedEdges/spline.C

@@ -36,7 +36,7 @@ Foam::spline::spline(const pointField& knotPoints)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::scalar Foam::spline::B(const scalar& tau) const
+Foam::scalar Foam::spline::B(const scalar tau) const
 {
     if (tau <= -2.0 || tau >= 2.0)
     {
@@ -60,7 +60,7 @@ Foam::scalar Foam::spline::B(const scalar& tau) const
     }
     else
     {
-        FatalErrorIn("spline::B(const scalar&)")
+        FatalErrorIn("spline::B(const scalar)")
             << "Programming error???, "
             << "tau = " << tau
             << abort(FatalError);

src/mesh/blockMesh/curvedEdges/spline.H

@@ -59,7 +59,7 @@ class spline
     // Private Member Functions
 
         //- Blending function for constructing spline
-        scalar B(const scalar&) const;
+        scalar B(const scalar) const;
 
         //- Disallow default bitwise copy construct
         spline(const spline&);

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

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)

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)