Merge commit 'OpenCFD/master' into olesenm

applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructor.H

@@ -151,18 +151,20 @@ public:
             const IOobject& fieldIoObject
         );
 
-        //- Reconstruct and write all volume fields
+        //- Reconstruct and write all/selected volume fields
         template<class Type>
         void reconstructFvVolumeFields
         (
-            const IOobjectList& objects
+            const IOobjectList& objects,
+            const HashSet<word>& selectedFields
         );
 
-        //- Reconstruct and write all volume fields
+        //- Reconstruct and write all/selected volume fields
         template<class Type>
         void reconstructFvSurfaceFields
         (
-            const IOobjectList& objects
+            const IOobjectList& objects,
+            const HashSet<word>& selectedFields
         );
 };
 

applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructorReconstructFields.C

@@ -131,7 +131,7 @@ Foam::fvFieldReconstructor::reconstructFvVolumeField
                 forAll(cp, faceI)
                 {
                     // Subtract one to take into account offsets for
-                    // face direction.  
+                    // face direction.
                     reverseAddressing[faceI] = cp[faceI] - 1 - curPatchStart;
                 }
 
@@ -151,7 +151,7 @@ Foam::fvFieldReconstructor::reconstructFvVolumeField
                 forAll(cp, faceI)
                 {
                     // Subtract one to take into account offsets for
-                    // face direction.  
+                    // face direction.
                     label curF = cp[faceI] - 1;
 
                     // Is the face on the boundary?
@@ -282,7 +282,7 @@ Foam::fvFieldReconstructor::reconstructFvSurfaceField
 
         // It is necessary to create a copy of the addressing array to
         // take care of the face direction offset trick.
-        // 
+        //
         {
             labelList curAddr(faceProcAddressing_[procI]);
 
@@ -342,7 +342,7 @@ Foam::fvFieldReconstructor::reconstructFvSurfaceField
                 forAll(cp, faceI)
                 {
                     // Subtract one to take into account offsets for
-                    // face direction.  
+                    // face direction.
                     reverseAddressing[faceI] = cp[faceI] - 1 - curPatchStart;
                 }
 
@@ -452,11 +452,12 @@ Foam::fvFieldReconstructor::reconstructFvSurfaceField
 }
 
 
-// Reconstruct and write all volume fields
+// Reconstruct and write all/selected volume fields
 template<class Type>
 void Foam::fvFieldReconstructor::reconstructFvVolumeFields
 (
-    const IOobjectList& objects
+    const IOobjectList& objects,
+    const HashSet<word>& selectedFields
 )
 {
     const word& fieldClassName =
@@ -468,27 +469,29 @@ void Foam::fvFieldReconstructor::reconstructFvVolumeFields
     {
         Info<< "    Reconstructing " << fieldClassName << "s\n" << endl;
 
-        for
-        (
-            IOobjectList::const_iterator fieldIter = fields.begin();
-            fieldIter != fields.end();
-            ++fieldIter
-        )
+        forAllConstIter(IOobjectList, fields, fieldIter)
         {
-            Info<< "        " << fieldIter()->name() << endl;
+            if
+            (
+                !selectedFields.size()
+             || selectedFields.found(fieldIter()->name())
+            )
+            {
+                Info<< "        " << fieldIter()->name() << endl;
 
-            reconstructFvVolumeField<Type>(*fieldIter())().write();
+                reconstructFvVolumeField<Type>(*fieldIter())().write();
+            }
         }
-
         Info<< endl;
     }
 }
 
-// Reconstruct and write all surface fields
+// Reconstruct and write all/selected surface fields
 template<class Type>
 void Foam::fvFieldReconstructor::reconstructFvSurfaceFields
 (
-    const IOobjectList& objects
+    const IOobjectList& objects,
+    const HashSet<word>& selectedFields
 )
 {
     const word& fieldClassName =
@@ -500,18 +503,19 @@ void Foam::fvFieldReconstructor::reconstructFvSurfaceFields
     {
         Info<< "    Reconstructing " << fieldClassName << "s\n" << endl;
 
-        for
-        (
-            IOobjectList::const_iterator fieldIter = fields.begin();
-            fieldIter != fields.end();
-            ++fieldIter
-        )
+        forAllConstIter(IOobjectList, fields, fieldIter)
         {
-            Info<< "        " << fieldIter()->name() << endl;
+            if
+            (
+                !selectedFields.size()
+             || selectedFields.found(fieldIter()->name())
+            )
+            {
+                Info<< "        " << fieldIter()->name() << endl;
 
-            reconstructFvSurfaceField<Type>(*fieldIter())().write();
+                reconstructFvSurfaceField<Type>(*fieldIter())().write();
+            }
         }
-
         Info<< endl;
     }
 }

applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C

@@ -48,9 +48,17 @@ int main(int argc, char *argv[])
     argList::noParallel();
     timeSelector::addOptions();
 #   include "addRegionOption.H"
+    argList::validOptions.insert("fields", "\"(list of fields)\"");
+
 #   include "setRootCase.H"
 #   include "createTime.H"
 
+    HashSet<word> selectedFields;
+    if (args.options().found("fields"))
+    {
+        IStringStream(args.options()["fields"])() >> selectedFields;
+    }
+
     // determine the processor count directly
     label nProcs = 0;
     while (dir(args.path()/(word("processor") + name(nProcs))))
@@ -184,13 +192,37 @@ int main(int argc, char *argv[])
                 procMeshes.boundaryProcAddressing()
             );
 
-            fvReconstructor.reconstructFvVolumeFields<scalar>(objects);
-            fvReconstructor.reconstructFvVolumeFields<vector>(objects);
-            fvReconstructor.reconstructFvVolumeFields<sphericalTensor>(objects);
-            fvReconstructor.reconstructFvVolumeFields<symmTensor>(objects);
-            fvReconstructor.reconstructFvVolumeFields<tensor>(objects);
+            fvReconstructor.reconstructFvVolumeFields<scalar>
+            (
+                objects,
+                selectedFields
+            );
+            fvReconstructor.reconstructFvVolumeFields<vector>
+            (
+                objects,
+                selectedFields
+            );
+            fvReconstructor.reconstructFvVolumeFields<sphericalTensor>
+            (
+                objects,
+                selectedFields
+            );
+            fvReconstructor.reconstructFvVolumeFields<symmTensor>
+            (
+                objects,
+                selectedFields
+            );
+            fvReconstructor.reconstructFvVolumeFields<tensor>
+            (
+                objects,
+                selectedFields
+            );
 
-            fvReconstructor.reconstructFvSurfaceFields<scalar>(objects);
+            fvReconstructor.reconstructFvSurfaceFields<scalar>
+            (
+                objects,
+                selectedFields
+            );
         }
         else
         {

applications/utilities/postProcessing/miscellaneous/execFlowFunctionObjects/execFlowFunctionObjects.C

@@ -139,7 +139,7 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
 
         IOobject LESPropertiesHeader
         (
-            "RASProperties",
+            "LESProperties",
             runTime.constant(),
             mesh,
             IOobject::MUST_READ,

src/ODE/ODE/ODE.H

@@ -34,7 +34,7 @@ Description
 #define ODE_H
 
 #include "scalarField.H"
-#include "Matrix.H"
+#include "scalarMatrices.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -80,7 +80,7 @@ public:
             const scalar x,
             const scalarField& y,
             scalarField& dfdx,
-            Matrix<scalar>& dfdy
+            scalarSquareMatrix& dfdy
         ) const = 0;
 };
 

src/ODE/ODESolvers/KRR4/KRR4.C

@@ -25,7 +25,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "KRR4.H"
-#include "simpleMatrix.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -36,11 +35,11 @@ namespace Foam
 {
     addToRunTimeSelectionTable(ODESolver, KRR4, ODE);
 
-const scalar 
+const scalar
     KRR4::safety = 0.9, KRR4::grow = 1.5, KRR4::pgrow = -0.25,
     KRR4::shrink = 0.5, KRR4::pshrink = (-1.0/3.0), KRR4::errcon = 0.1296;
 
-const scalar 
+const scalar
     KRR4::gamma = 1.0/2.0,
     KRR4::a21 = 2.0, KRR4::a31 = 48.0/25.0, KRR4::a32 = 6.0/25.0,
     KRR4::c21 = -8.0, KRR4::c31 = 372.0/25.0, KRR4::c32 = 12.0/5.0,
@@ -81,8 +80,8 @@ void Foam::KRR4::solve
     const ODE& ode,
     scalar& x,
     scalarField& y,
-    scalarField& dydx, 
-    const scalar eps, 
+    scalarField& dydx,
+    const scalar eps,
     const scalarField& yScale,
     const scalar hTry,
     scalar& hDid,
@@ -109,14 +108,14 @@ void Foam::KRR4::solve
             a_[i][i] += 1.0/(gamma*h);
         }
 
-        simpleMatrix<scalar>::LUDecompose(a_, pivotIndices_);
+        LUDecompose(a_, pivotIndices_);
 
         for (register label i=0; i<n_; i++)
         {
             g1_[i] = dydxTemp_[i] + h*c1X*dfdx_[i];
         }
 
-        simpleMatrix<scalar>::LUBacksubstitute(a_, pivotIndices_, g1_);
+        LUBacksubstitute(a_, pivotIndices_, g1_);
 
         for (register label i=0; i<n_; i++)
         {
@@ -131,7 +130,7 @@ void Foam::KRR4::solve
             g2_[i] = dydx_[i] + h*c2X*dfdx_[i] + c21*g1_[i]/h;
         }
 
-        simpleMatrix<scalar>::LUBacksubstitute(a_, pivotIndices_, g2_);
+        LUBacksubstitute(a_, pivotIndices_, g2_);
 
         for (register label i=0; i<n_; i++)
         {
@@ -146,7 +145,7 @@ void Foam::KRR4::solve
             g3_[i] = dydx[i] + h*c3X*dfdx_[i] + (c31*g1_[i] + c32*g2_[i])/h;
         }
 
-        simpleMatrix<scalar>::LUBacksubstitute(a_, pivotIndices_, g3_);
+        LUBacksubstitute(a_, pivotIndices_, g3_);
 
         for (register label i=0; i<n_; i++)
         {
@@ -154,7 +153,7 @@ void Foam::KRR4::solve
                 + (c41*g1_[i] + c42*g2_[i] + c43*g3_[i])/h;
         }
 
-        simpleMatrix<scalar>::LUBacksubstitute(a_, pivotIndices_, g4_);
+        LUBacksubstitute(a_, pivotIndices_, g4_);
 
         for (register label i=0; i<n_; i++)
         {

src/ODE/ODESolvers/KRR4/KRR4.H

@@ -62,15 +62,15 @@ class KRR4
         mutable scalarField g4_;
         mutable scalarField yErr_;
         mutable scalarField dfdx_;
-        mutable Matrix<scalar> dfdy_;
-        mutable Matrix<scalar> a_;
+        mutable scalarSquareMatrix dfdy_;
+        mutable scalarSquareMatrix a_;
         mutable labelList pivotIndices_;
 
         static const int maxtry = 40;
 
         static const scalar safety, grow, pgrow, shrink, pshrink, errcon;
 
-        static const scalar 
+        static const scalar
             gamma,
             a21, a31, a32,
             c21, c31, c32, c41, c42, c43,

src/ODE/ODESolvers/SIBS/SIBS.H

@@ -60,8 +60,8 @@ class SIBS
         static const scalar safe1, safe2, redMax, redMin, scaleMX;
 
         mutable scalarField a_;
-        mutable Matrix<scalar> alpha_;
-        mutable Matrix<scalar> d_p_;
+        mutable scalarSquareMatrix alpha_;
+        mutable scalarSquareMatrix d_p_;
         mutable scalarField x_p_;
         mutable scalarField err_;
 
@@ -69,7 +69,7 @@ class SIBS
         mutable scalarField ySeq_;
         mutable scalarField yErr_;
         mutable scalarField dfdx_;
-        mutable Matrix<scalar> dfdy_;
+        mutable scalarSquareMatrix dfdy_;
 
         mutable label first_, kMax_, kOpt_;
         mutable scalar epsOld_, xNew_;
@@ -82,7 +82,7 @@ void SIMPR
     const scalarField& y,
     const scalarField& dydx,
     const scalarField& dfdx,
-    const Matrix<scalar>& dfdy, 
+    const scalarSquareMatrix& dfdy,
     const scalar deltaX,
     const label nSteps,
     scalarField& yEnd
@@ -97,7 +97,7 @@ void polyExtrapolate
     scalarField& yz,
     scalarField& dy,
     scalarField& x_p,
-    Matrix<scalar>& d_p
+    scalarSquareMatrix& d_p
 ) const;
 
 

src/ODE/ODESolvers/SIBS/SIMPR.C

@@ -25,7 +25,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "SIBS.H"
-#include "simpleMatrix.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -36,7 +35,7 @@ void Foam::SIBS::SIMPR
     const scalarField& y,
     const scalarField& dydx,
     const scalarField& dfdx,
-    const Matrix<scalar>& dfdy,
+    const scalarSquareMatrix& dfdy,
     const scalar deltaX,
     const label nSteps,
     scalarField& yEnd
@@ -44,7 +43,7 @@ void Foam::SIBS::SIMPR
 {
     scalar h = deltaX/nSteps;
 
-    Matrix<scalar> a(n_, n_);
+    scalarSquareMatrix a(n_);
     for (register label i=0; i<n_; i++)
     {
         for (register label j=0; j<n_; j++)
@@ -55,14 +54,14 @@ void Foam::SIBS::SIMPR
     }
 
     labelList pivotIndices(n_);
-    simpleMatrix<scalar>::LUDecompose(a, pivotIndices);
+    LUDecompose(a, pivotIndices);
 
     for (register label i=0; i<n_; i++)
     {
         yEnd[i] = h*(dydx[i] + h*dfdx[i]);
     }
 
-    simpleMatrix<scalar>::LUBacksubstitute(a, pivotIndices, yEnd);
+    LUBacksubstitute(a, pivotIndices, yEnd);
 
     scalarField del(yEnd);
     scalarField ytemp(n_);
@@ -83,7 +82,7 @@ void Foam::SIBS::SIMPR
             yEnd[i] = h*yEnd[i] - del[i];
         }
 
-        simpleMatrix<scalar>::LUBacksubstitute(a, pivotIndices, yEnd);
+        LUBacksubstitute(a, pivotIndices, yEnd);
 
         for (register label i=0; i<n_; i++)
         {
@@ -99,7 +98,7 @@ void Foam::SIBS::SIMPR
         yEnd[i] = h*yEnd[i] - del[i];
     }
 
-    simpleMatrix<scalar>::LUBacksubstitute(a, pivotIndices, yEnd);
+    LUBacksubstitute(a, pivotIndices, yEnd);
 
     for (register label i=0; i<n_; i++)
     {

src/ODE/ODESolvers/SIBS/polyExtrapolate.C

@@ -36,7 +36,7 @@ void Foam::SIBS::polyExtrapolate
     scalarField& yz,
     scalarField& dy,
     scalarField& x,
-    Matrix<scalar>& d
+    scalarSquareMatrix& d
 ) const
 {
     label n = yz.size();

src/OpenFOAM/Make/files

@@ -177,8 +177,9 @@ dimensionedTypes/dimensionedTensor/dimensionedTensor.C
 
 matrices/solution/solution.C
 
-scalarMatrix = matrices/scalarMatrix
-$(scalarMatrix)/scalarMatrix.C
+scalarMatrices = matrices/scalarMatrices
+$(scalarMatrices)/scalarMatrices.C
+$(scalarMatrices)/SVD/SVD.C
 
 LUscalarMatrix = matrices/LUscalarMatrix
 $(LUscalarMatrix)/LUscalarMatrix.C

src/OpenFOAM/db/IOstreams/Pstreams/PstreamCommsStruct.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 "Pstream.H"

src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C

@@ -0,0 +1,100 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "DiagonalMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+template<class Form>
+Foam::DiagonalMatrix<Type>::DiagonalMatrix(const Matrix<Form, Type>& a)
+:
+    List<Type>(min(a.n(), a.m()))
+{
+    forAll(*this, i)
+    {
+        this->operator[](i) = a[i][i];
+    }
+}
+
+
+template<class Type>
+Foam::DiagonalMatrix<Type>::DiagonalMatrix(const label size)
+:
+    List<Type>(size)
+{}
+
+
+template<class Type>
+Foam::DiagonalMatrix<Type>::DiagonalMatrix(const label size, const Type& val)
+:
+    List<Type>(size, val)
+{}
+
+
+template<class Type>
+Foam::DiagonalMatrix<Type>& Foam::DiagonalMatrix<Type>::invert()
+{
+    forAll(*this, i)
+    {
+        Type x = this->operator[](i);
+        if (mag(x) < VSMALL)
+        {
+            this->operator[](i) = Type(0);
+        }
+        else
+        {
+            this->operator[](i) = Type(1)/x;
+        }
+    }
+
+    return this;
+}
+
+
+template<class Type>
+Foam::DiagonalMatrix<Type> Foam::inv(const DiagonalMatrix<Type>& A)
+{
+    DiagonalMatrix<Type> Ainv = A;
+
+    forAll(A, i)
+    {
+        Type x = A[i];
+        if (mag(x) < VSMALL)
+        {
+            Ainv[i] = Type(0);
+        }
+        else
+        {
+            Ainv[i] = Type(1)/x;
+        }
+    }
+
+    return Ainv;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H

@@ -0,0 +1,103 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::DiagonalMatrix<Type>
+
+Description
+    DiagonalMatrix<Type> is a 2D diagonal matrix of objects
+    of type Type, size nxn
+
+SourceFiles
+    DiagonalMatrix.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef DiagonalMatrix_H
+#define DiagonalMatrix_H
+
+#include "List.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * *  * * * * * * Class Forward declaration  * * * * * * * * * * * //
+
+template<class Form, class Type> class Matrix;
+
+/*---------------------------------------------------------------------------*\
+                           Class DiagonalMatrix Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class DiagonalMatrix
+:
+    public List<Type>
+{
+public:
+
+    // Constructors
+
+        //- Construct from diagonal component of a Matrix
+        template<class Form>
+        DiagonalMatrix<Type>(const Matrix<Form, Type>&);
+
+        //- Construct empty from size
+        DiagonalMatrix<Type>(const label size);
+
+        //- Construct from size and a value
+        DiagonalMatrix<Type>(const label, const Type&);
+
+
+    // Member functions
+
+        //- Invert the diaganol matrix and return itself
+        DiagonalMatrix<Type>& invert();
+};
+
+
+// Global functions
+
+//- Return the diagonal Matrix inverse
+template<class Type>
+DiagonalMatrix<Type> inv(const DiagonalMatrix<Type>&);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "DiagonalMatrix.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C

@@ -31,9 +31,9 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::LUscalarMatrix::LUscalarMatrix(const Matrix<scalar>& matrix)
+Foam::LUscalarMatrix::LUscalarMatrix(const scalarSquareMatrix& matrix)
 :
-    scalarMatrix(matrix),
+    scalarSquareMatrix(matrix),
     pivotIndices_(n())
 {
     LUDecompose(*this, pivotIndices_);
@@ -101,7 +101,7 @@ Foam::LUscalarMatrix::LUscalarMatrix
                 nCells += lduMatrices[i].size();
             }
 
-            Matrix<scalar> m(nCells, nCells, 0.0);
+            scalarSquareMatrix m(nCells, nCells, 0.0);
             transfer(m);
             convert(lduMatrices);
         }
@@ -109,7 +109,7 @@ Foam::LUscalarMatrix::LUscalarMatrix
     else
     {
         label nCells = ldum.lduAddr().size();
-        Matrix<scalar> m(nCells, nCells, 0.0);
+        scalarSquareMatrix m(nCells, nCells, 0.0);
         transfer(m);
         convert(ldum, interfaceCoeffs, interfaces);
     }

src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H

@@ -36,7 +36,7 @@ SourceFiles
 #ifndef LUscalarMatrix_H
 #define LUscalarMatrix_H
 
-#include "scalarMatrix.H"
+#include "scalarMatrices.H"
 #include "labelList.H"
 #include "FieldField.H"
 #include "lduInterfaceFieldPtrsList.H"
@@ -55,7 +55,7 @@ class procLduMatrix;
 
 class LUscalarMatrix
 :
-    public scalarMatrix
+    public scalarSquareMatrix
 {
     // Private data
 
@@ -78,7 +78,7 @@ class LUscalarMatrix
         void convert(const PtrList<procLduMatrix>& lduMatrices);
 
 
-        //- Print the ratio of the mag-sum of the off-diagonal coefficients 
+        //- Print the ratio of the mag-sum of the off-diagonal coefficients
         //  to the mag-diagonal
         void printDiagonalDominance() const;
 
@@ -87,8 +87,8 @@ public:
 
     // Constructors
 
-        //- Construct from Matrix<scalar> and perform LU decomposition
-        LUscalarMatrix(const Matrix<scalar>&);
+        //- Construct from scalarSquareMatrix and perform LU decomposition
+        LUscalarMatrix(const scalarSquareMatrix&);
 
         //- Construct from lduMatrix and perform LU decomposition
         LUscalarMatrix

src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C

@@ -28,16 +28,16 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class T>
-void Foam::LUscalarMatrix::solve(Field<T>& sourceSol) const
+template<class Type>
+void Foam::LUscalarMatrix::solve(Field<Type>& sourceSol) const
 {
     if (Pstream::parRun())
     {
-        Field<T> completeSourceSol(n());
+        Field<Type> completeSourceSol(n());
 
         if (Pstream::master())
         {
-            typename Field<T>::subField
+            typename Field<Type>::subField
             (
                 completeSourceSol,
                 sourceSol.size()
@@ -58,7 +58,7 @@ void Foam::LUscalarMatrix::solve(Field<T>& sourceSol) const
                     (
                         &(completeSourceSol[procOffsets_[slave]])
                     ),
-                    (procOffsets_[slave + 1] - procOffsets_[slave])*sizeof(T)
+                    (procOffsets_[slave + 1] - procOffsets_[slave])*sizeof(Type)
                 );
             }
         }
@@ -77,7 +77,7 @@ void Foam::LUscalarMatrix::solve(Field<T>& sourceSol) const
         {
             LUBacksubstitute(*this, pivotIndices_, completeSourceSol);
 
-            sourceSol = typename Field<T>::subField
+            sourceSol = typename Field<Type>::subField
             (
                 completeSourceSol,
                 sourceSol.size()
@@ -98,7 +98,7 @@ void Foam::LUscalarMatrix::solve(Field<T>& sourceSol) const
                     (
                         &(completeSourceSol[procOffsets_[slave]])
                     ),
-                    (procOffsets_[slave + 1] - procOffsets_[slave])*sizeof(T)
+                    (procOffsets_[slave + 1] - procOffsets_[slave])*sizeof(Type)
                 );
             }
         }

src/OpenFOAM/matrices/Matrix/Matrix.C

@@ -28,13 +28,13 @@ License
 
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
-template<class T>
-void Foam::Matrix<T>::allocate()
+template<class Form, class Type>
+void Foam::Matrix<Form, Type>::allocate()
 {
     if (n_ && m_)
     {
-        v_ = new T*[n_];
-        v_[0] = new T[n_*m_];
+        v_ = new Type*[n_];
+        v_[0] = new Type[n_*m_];
 
         for (register label i=1; i<n_; i++)
         {
@@ -46,12 +46,12 @@ void Foam::Matrix<T>::allocate()
 
 // * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * * //
 
-template<class T>
-Foam::Matrix<T>::~Matrix()
+template<class Form, class Type>
+Foam::Matrix<Form, Type>::~Matrix()
 {
     if (v_)
     {
-		delete[] (v_[0]);
+        delete[] (v_[0]);
         delete[] v_;
     }
 }
@@ -59,16 +59,16 @@ Foam::Matrix<T>::~Matrix()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class T>
-const Foam::Matrix<T>& Foam::Matrix<T>::null()
+template<class Form, class Type>
+const Foam::Matrix<Form, Type>& Foam::Matrix<Form, Type>::null()
 {
-    Matrix<T>* nullPtr = reinterpret_cast<Matrix<T>*>(NULL);
+    Matrix<Form, Type>* nullPtr = reinterpret_cast<Matrix<Form, Type>*>(NULL);
     return *nullPtr;
 }
 
 
-template<class T>
-Foam::Matrix<T>::Matrix(const label n, const label m)
+template<class Form, class Type>
+Foam::Matrix<Form, Type>::Matrix(const label n, const label m)
 :
     n_(n),
     m_(m),
@@ -76,7 +76,7 @@ Foam::Matrix<T>::Matrix(const label n, const label m)
 {
     if (n_ < 0 || m_ < 0)
     {
-        FatalErrorIn("Matrix<T>::Matrix(const label n, const label m)")
+        FatalErrorIn("Matrix<Form, Type>::Matrix(const label n, const label m)")
             << "bad n, m " << n_ << ", " << m_
             << abort(FatalError);
     }
@@ -85,8 +85,8 @@ Foam::Matrix<T>::Matrix(const label n, const label m)
 }
 
 
-template<class T>
-Foam::Matrix<T>::Matrix(const label n, const label m, const T& a)
+template<class Form, class Type>
+Foam::Matrix<Form, Type>::Matrix(const label n, const label m, const Type& a)
 :
     n_(n),
     m_(m),
@@ -96,7 +96,7 @@ Foam::Matrix<T>::Matrix(const label n, const label m, const T& a)
     {
         FatalErrorIn
         (
-            "Matrix<T>::Matrix(const label n, const label m, const T&)"
+            "Matrix<Form, Type>::Matrix(const label n, const label m, const T&)"
         )   << "bad n, m " << n_ << ", " << m_
             << abort(FatalError);
     }
@@ -105,7 +105,7 @@ Foam::Matrix<T>::Matrix(const label n, const label m, const T& a)
 
     if (v_)
     {
-        T* v = v_[0];
+        Type* v = v_[0];
 
         label nm = n_*m_;
 
@@ -117,8 +117,8 @@ Foam::Matrix<T>::Matrix(const label n, const label m, const T& a)
 }
 
 
-template<class T>
-Foam::Matrix<T>::Matrix(const Matrix<T>& a)
+template<class Form, class Type>
+Foam::Matrix<Form, Type>::Matrix(const Matrix<Form, Type>& a)
 :
     n_(a.n_),
     m_(a.m_),
@@ -127,8 +127,8 @@ Foam::Matrix<T>::Matrix(const Matrix<T>& a)
     if (a.v_)
     {
         allocate();
-        T* v = v_[0];
-        const T* av = a.v_[0];
+        Type* v = v_[0];
+        const Type* av = a.v_[0];
 
         label nm = n_*m_;
         for (register label i=0; i<nm; i++)
@@ -138,13 +138,13 @@ Foam::Matrix<T>::Matrix(const Matrix<T>& a)
     }
 }
 
-        
-template<class T>
-void Foam::Matrix<T>::clear()
+
+template<class Form, class Type>
+void Foam::Matrix<Form, Type>::clear()
 {
     if (v_)
     {
-		delete[] (v_[0]);
+        delete[] (v_[0]);
         delete[] v_;
     }
     n_ = 0;
@@ -153,8 +153,8 @@ void Foam::Matrix<T>::clear()
 }
 
 
-template<class T>
-void Foam::Matrix<T>::transfer(Matrix<T>& a)
+template<class Form, class Type>
+void Foam::Matrix<Form, Type>::transfer(Matrix<Form, Type>& a)
 {
     clear();
 
@@ -169,14 +169,32 @@ void Foam::Matrix<T>::transfer(Matrix<T>& a)
 }
 
 
+template<class Form, class Type>
+Form Foam::Matrix<Form, Type>::T() const
+{
+    const Matrix<Form, Type>& A = *this;
+    Form At(m(), n());
+
+    for (register label i=0; i<n(); i++)
+    {
+        for (register label j=0; j<m(); j++)
+        {
+            At[j][i] = A[i][j];
+        }
+    }
+
+    return At;
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
-template<class T>
-void Foam::Matrix<T>::operator=(const T& t)
+template<class Form, class Type>
+void Foam::Matrix<Form, Type>::operator=(const Type& t)
 {
     if (v_)
     {
-        T* v = v_[0];
+        Type* v = v_[0];
 
         label nm = n_*m_;
         for (register label i=0; i<nm; i++)
@@ -188,12 +206,12 @@ void Foam::Matrix<T>::operator=(const T& t)
 
 
 // Assignment operator. Takes linear time.
-template<class T>
-void Foam::Matrix<T>::operator=(const Matrix<T>& a)
+template<class Form, class Type>
+void Foam::Matrix<Form, Type>::operator=(const Matrix<Form, Type>& a)
 {
     if (this == &a)
     {
-        FatalErrorIn("Matrix<T>::operator=(const Matrix<T>&)")
+        FatalErrorIn("Matrix<Form, Type>::operator=(const Matrix<Form, Type>&)")
             << "attempted assignment to self"
             << abort(FatalError);
     }
@@ -204,12 +222,12 @@ void Foam::Matrix<T>::operator=(const Matrix<T>& a)
         n_ = a.n_;
         m_ = a.m_;
         allocate();
-    } 
+    }
 
     if (v_)
     {
-        T* v = v_[0];
-        const T* av = a.v_[0];
+        Type* v = v_[0];
+        const Type* av = a.v_[0];
 
         label nm = n_*m_;
         for (register label i=0; i<nm; i++)
@@ -222,15 +240,15 @@ void Foam::Matrix<T>::operator=(const Matrix<T>& a)
 
 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 
-template<class T>
-const T& Foam::max(const Matrix<T>& a)
+template<class Form, class Type>
+const Type& Foam::max(const Matrix<Form, Type>& a)
 {
     label nm = a.n_*a.m_;
 
     if (nm)
     {
         label curMaxI = 0;
-        const T* v = a.v_[0];
+        const Type* v = a.v_[0];
 
         for (register label i=1; i<nm; i++)
         {
@@ -244,7 +262,7 @@ const T& Foam::max(const Matrix<T>& a)
     }
     else
     {
-        FatalErrorIn("max(const Matrix<T>&)")
+        FatalErrorIn("max(const Matrix<Form, Type>&)")
             << "matrix is empty"
             << abort(FatalError);
 
@@ -254,15 +272,15 @@ const T& Foam::max(const Matrix<T>& a)
 }
 
 
-template<class T>
-const T& Foam::min(const Matrix<T>& a)
+template<class Form, class Type>
+const Type& Foam::min(const Matrix<Form, Type>& a)
 {
     label nm = a.n_*a.m_;
 
     if (nm)
     {
         label curMinI = 0;
-        const T* v = a.v_[0];
+        const Type* v = a.v_[0];
 
         for (register label i=1; i<nm; i++)
         {
@@ -276,7 +294,7 @@ const T& Foam::min(const Matrix<T>& a)
     }
     else
     {
-        FatalErrorIn("min(const Matrix<T>&)")
+        FatalErrorIn("min(const Matrix<Form, Type>&)")
             << "matrix is empty"
             << abort(FatalError);
 
@@ -288,15 +306,15 @@ const T& Foam::min(const Matrix<T>& a)
 
 // * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
 
-template<class T>
-Foam::Matrix<T> Foam::operator-(const Matrix<T>& a)
+template<class Form, class Type>
+Form Foam::operator-(const Matrix<Form, Type>& a)
 {
-    Matrix<T> na(a.n_, a.m_);
+    Form na(a.n_, a.m_);
 
     if (a.n_ && a.m_)
     {
-        T* nav = na.v_[0];
-        const T* av = a.v_[0];
+        Type* nav = na.v_[0];
+        const Type* av = a.v_[0];
 
         label nm = a.n_*a.m_;
         for (register label i=0; i<nm; i++)
@@ -309,30 +327,34 @@ Foam::Matrix<T> Foam::operator-(const Matrix<T>& a)
 }
 
 
-template<class T>
-Foam::Matrix<T> Foam::operator+(const Matrix<T>& a, const Matrix<T>& b)
+template<class Form, class Type>
+Form Foam::operator+(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
 {
     if (a.n_ != b.n_)
     {
-        FatalErrorIn("Matrix<T>::operator+(const Matrix<T>&, const Matrix<T>&)")
-            << "attempted add matrices with different number of rows: "
+        FatalErrorIn
+        (
+            "Matrix<Form, Type>::operator+(const Matrix<Form, Type>&, const Matrix<Form, Type>&)"
+        )   << "attempted add matrices with different number of rows: "
             << a.n_ << ", " << b.n_
             << abort(FatalError);
     }
 
     if (a.m_ != b.m_)
     {
-        FatalErrorIn("Matrix<T>::operator+(const Matrix<T>&, const Matrix<T>&)")
-            << "attempted add matrices with different number of columns: "
+        FatalErrorIn
+        (
+            "Matrix<Form, Type>::operator+(const Matrix<Form, Type>&, const Matrix<Form, Type>&)"
+        )   << "attempted add matrices with different number of columns: "
             << a.m_ << ", " << b.m_
             << abort(FatalError);
     }
 
-    Matrix<T> ab(a.n_, a.m_);
+    Form ab(a.n_, a.m_);
 
-    T* abv = ab.v_[0];
-    const T* av = a.v_[0];
-    const T* bv = b.v_[0];
+    Type* abv = ab.v_[0];
+    const Type* av = a.v_[0];
+    const Type* bv = b.v_[0];
 
     label nm = a.n_*a.m_;;
     for (register label i=0; i<nm; i++)
@@ -344,30 +366,34 @@ Foam::Matrix<T> Foam::operator+(const Matrix<T>& a, const Matrix<T>& b)
 }
 
 
-template<class T>
-Foam::Matrix<T> Foam::operator-(const Matrix<T>& a, const Matrix<T>& b)
+template<class Form, class Type>
+Form Foam::operator-(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
 {
     if (a.n_ != b.n_)
     {
-        FatalErrorIn("Matrix<T>::operator-(const Matrix<T>&, const Matrix<T>&)")
-            << "attempted add matrices with different number of rows: "
+        FatalErrorIn
+        (
+            "Matrix<Form, Type>::operator-(const Matrix<Form, Type>&, const Matrix<Form, Type>&)"
+        )   << "attempted add matrices with different number of rows: "
             << a.n_ << ", " << b.n_
             << abort(FatalError);
     }
 
     if (a.m_ != b.m_)
     {
-        FatalErrorIn("Matrix<T>::operator-(const Matrix<T>&, const Matrix<T>&)")
-            << "attempted add matrices with different number of columns: "
+        FatalErrorIn
+        (
+            "Matrix<Form, Type>::operator-(const Matrix<Form, Type>&, const Matrix<Form, Type>&)"
+        )   << "attempted add matrices with different number of columns: "
             << a.m_ << ", " << b.m_
             << abort(FatalError);
     }
 
-    Matrix<T> ab(a.n_, a.m_);
+    Form ab(a.n_, a.m_);
 
-    T* abv = ab.v_[0];
-    const T* av = a.v_[0];
-    const T* bv = b.v_[0];
+    Type* abv = ab.v_[0];
+    const Type* av = a.v_[0];
+    const Type* bv = b.v_[0];
 
     label nm = a.n_*a.m_;;
     for (register label i=0; i<nm; i++)
@@ -379,17 +405,17 @@ Foam::Matrix<T> Foam::operator-(const Matrix<T>& a, const Matrix<T>& b)
 }
 
 
-template<class T>
-Foam::Matrix<T> Foam::operator*(const scalar s, const Matrix<T>& a)
+template<class Form, class Type>
+Form Foam::operator*(const scalar s, const Matrix<Form, Type>& a)
 {
-    Matrix<T> sa(a.n_, a.m_);
+    Form sa(a.n_, a.m_);
 
     if (a.n_ && a.m_)
     {
-        T* sav = sa.v_[0];
-        const T* av = a.v_[0];
+        Type* sav = sa.v_[0];
+        const Type* av = a.v_[0];
 
-        label nm = a.n_*a.m_;;
+        label nm = a.n_*a.m_;
         for (register label i=0; i<nm; i++)
         {
             sav[i] = s*av[i];

src/OpenFOAM/matrices/Matrix/Matrix.H

@@ -51,25 +51,26 @@ namespace Foam
 
 // Forward declaration of friend functions and operators
 
-template<class T> class Matrix;
+template<class Form, class Type> class Matrix;
 
-template<class T> const T& max(const Matrix<T>&);
-template<class T> const T& min(const Matrix<T>&);
+template<class Form, class Type> Istream& operator>>
+(
+    Istream&,
+    Matrix<Form, Type>&
+);
 
-template<class T> Matrix<T> operator-(const Matrix<T>&);
-template<class T> Matrix<T> operator+(const Matrix<T>&, const Matrix<T>&);
-template<class T> Matrix<T> operator-(const Matrix<T>&, const Matrix<T>&);
-template<class T> Matrix<T> operator*(const scalar, const Matrix<T>&);
-
-template<class T> Istream& operator>>(Istream&, Matrix<T>&);
-template<class T> Ostream& operator<<(Ostream&, const Matrix<T>&);
+template<class Form, class Type> Ostream& operator<<
+(
+    Ostream&,
+    const Matrix<Form, Type>&
+);
 
 
 /*---------------------------------------------------------------------------*\
                            Class Matrix Declaration
 \*---------------------------------------------------------------------------*/
 
-template<class T>
+template<class Form, class Type>
 class Matrix
 {
     // Private data
@@ -78,7 +79,7 @@ class Matrix
         label n_, m_;
 
         //- Row pointers
-        T** __restrict__ v_;
+        Type** __restrict__ v_;
 
         //- Allocate the storage for the row-pointers and the data
         //  and set the row pointers
@@ -97,16 +98,16 @@ public:
 
         //- Construct with given number of rows and columns
         //  and value for all elements.
-        Matrix(const label n, const label m, const T&);
+        Matrix(const label n, const label m, const Type&);
 
         //- Copy constructor.
-        Matrix(const Matrix<T>&);
+        Matrix(const Matrix<Form, Type>&);
 
         //- Construct from Istream.
         Matrix(Istream&);
 
         //- Clone
-        inline autoPtr<Matrix<T> > clone() const;
+        inline autoPtr<Matrix<Form, Type> > clone() const;
 
 
     // Destructor
@@ -117,7 +118,7 @@ public:
     // Member functions
 
         //- Return a null Matrix
-        static const Matrix<T>& null();
+        static const Matrix<Form, Type>& null();
 
 
         // Access
@@ -148,48 +149,64 @@ public:
 
             //- Transfer the contents of the argument Matrix into this Matrix
             //  and annull the argument Matrix.
-            void transfer(Matrix<T>&);
+            void transfer(Matrix<Form, Type>&);
+
+
+        //- Return the transpose of the matrix
+        Form T() const;
 
 
     // Member operators
 
         //- Return subscript-checked element of Matrix.
-        inline T* operator[](const label);
+        inline Type* operator[](const label);
 
         //- Return subscript-checked element of constant Matrix.
-        inline const T* operator[](const label) const;
+        inline const Type* operator[](const label) const;
 
         //- Assignment operator. Takes linear time.
-        void operator=(const Matrix<T>&);
+        void operator=(const Matrix<Form, Type>&);
 
         //- Assignment of all entries to the given value
-        void operator=(const T&);
-
-
-    // Friend functions
-
-        friend const T& max<T>(const Matrix<T>&);
-        friend const T& min<T>(const Matrix<T>&);
-
-
-    // Friend operators
-
-        friend Matrix<T> operator-<T>(const Matrix<T>&);
-        friend Matrix<T> operator+<T>(const Matrix<T>&, const Matrix<T>&);
-        friend Matrix<T> operator-<T>(const Matrix<T>&, const Matrix<T>&);
-        friend Matrix<T> operator*<T>(const scalar, const Matrix<T>&);
+        void operator=(const Type&);
 
 
     // IOstream operators
 
         //- Read Matrix from Istream, discarding contents of existing Matrix.
-        friend Istream& operator>> <T>(Istream&, Matrix<T>&);
+        friend Istream& operator>> <Form, Type>(Istream&, Matrix<Form, Type>&);
 
         // Write Matrix to Ostream.
-        friend Ostream& operator<< <T>(Ostream&, const Matrix<T>&);
+        friend Ostream& operator<< <Form, Type>(Ostream&, const Matrix<Form, Type>&);
 };
 
 
+// Global functions and operators
+
+template<class Form, class Type> const Type& max(const Matrix<Form, Type>&);
+template<class Form, class Type> const Type& min(const Matrix<Form, Type>&);
+
+template<class Form, class Type> Form operator-(const Matrix<Form, Type>&);
+
+template<class Form, class Type> Form operator+
+(
+    const Matrix<Form, Type>&,
+    const Matrix<Form, Type>&
+);
+
+template<class Form, class Type> Form operator-
+(
+    const Matrix<Form, Type>&,
+    const Matrix<Form, Type>&
+);
+
+template<class Form, class Type> Form operator*
+(
+    const scalar,
+    const Matrix<Form, Type>&
+);
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/matrices/Matrix/MatrixI.H

@@ -26,8 +26,8 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-template<class T>
-inline Foam::Matrix<T>::Matrix()
+template<class Form, class Type>
+inline Foam::Matrix<Form, Type>::Matrix()
 :
     n_(0),
     m_(0),
@@ -35,71 +35,67 @@ inline Foam::Matrix<T>::Matrix()
 {}
 
 
-template<class T>
-inline Foam::autoPtr<Foam::Matrix<T> > Foam::Matrix<T>::clone() const
+template<class Form, class Type>
+inline Foam::autoPtr<Foam::Matrix<Form, Type> > Foam::Matrix<Form, Type>::clone() const
 {
-    return autoPtr<Matrix<T> >(new Matrix<T>(*this));
+    return autoPtr<Matrix<Form, Type> >(new Matrix<Form, Type>(*this));
 }
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 //- Return the number of rows
-template<class T>
-inline Foam::label Foam::Matrix<T>::n() const
+template<class Form, class Type>
+inline Foam::label Foam::Matrix<Form, Type>::n() const
 {
     return n_;
 }
 
 
-//- Return the number of columns
-template<class T>
-inline Foam::label Foam::Matrix<T>::m() const
+template<class Form, class Type>
+inline Foam::label Foam::Matrix<Form, Type>::m() const
 {
     return m_;
 }
 
 
-//- Return the number of columns
-template<class T>
-inline Foam::label Foam::Matrix<T>::size() const
+template<class Form, class Type>
+inline Foam::label Foam::Matrix<Form, Type>::size() const
 {
     return n_*m_;
 }
 
 
-// Check index i is within valid range (0 ... n-1).
-template<class T>
-inline void Foam::Matrix<T>::checki(const label i) const
+template<class Form, class Type>
+inline void Foam::Matrix<Form, Type>::checki(const label i) const
 {
     if (!n_)
     {
-        FatalErrorIn("Matrix<T>::checki(const label)")
+        FatalErrorIn("Matrix<Form, Type>::checki(const label)")
             << "attempt to access element from zero sized row"
             << abort(FatalError);
     }
     else if (i<0 || i>=n_)
     {
-        FatalErrorIn("Matrix<T>::checki(const label)")
+        FatalErrorIn("Matrix<Form, Type>::checki(const label)")
             << "index " << i << " out of range 0 ... " << n_-1
             << abort(FatalError);
     }
 }
 
 
-// Check index j is within valid range (0 ... n-1).
-template<class T>
-inline void Foam::Matrix<T>::checkj(const label j) const
+template<class Form, class Type>
+inline void Foam::Matrix<Form, Type>::checkj(const label j) const
 {
     if (!m_)
     {
-        FatalErrorIn("Matrix<T>::checkj(const label)")
+        FatalErrorIn("Matrix<Form, Type>::checkj(const label)")
             << "attempt to access element from zero sized column"
             << abort(FatalError);
     }
     else if (j<0 || j>=m_)
     {
-        FatalErrorIn("Matrix<T>::checkj(const label)")
+        FatalErrorIn("Matrix<Form, Type>::checkj(const label)")
             << "index " << j << " out of range 0 ... " << m_-1
             << abort(FatalError);
     }
@@ -108,9 +104,8 @@ inline void Foam::Matrix<T>::checkj(const label j) const
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
-// Return subscript-checked element access
-template<class T>
-inline T* Foam::Matrix<T>::operator[](const label i)
+template<class Form, class Type>
+inline Type* Foam::Matrix<Form, Type>::operator[](const label i)
 {
 #   ifdef FULLDEBUG
     checki(i);
@@ -119,9 +114,8 @@ inline T* Foam::Matrix<T>::operator[](const label i)
 }
 
 
-// Return subscript-checked const element access
-template<class T>
-inline const T* Foam::Matrix<T>::operator[](const label i) const
+template<class Form, class Type>
+inline const Type* Foam::Matrix<Form, Type>::operator[](const label i) const
 {
 #   ifdef FULLDEBUG
     checki(i);

src/OpenFOAM/matrices/Matrix/MatrixIO.C

@@ -32,8 +32,8 @@ License
 
 // * * * * * * * * * * * * * * * Ostream Operator *  * * * * * * * * * * * * //
 
-template<class T>
-Foam::Matrix<T>::Matrix(Istream& is)
+template<class Form, class Type>
+Foam::Matrix<Form, Type>::Matrix(Istream& is)
 :
     n_(0),
     m_(0),
@@ -43,17 +43,17 @@ Foam::Matrix<T>::Matrix(Istream& is)
 }
 
 
-template<class T>
-Foam::Istream& Foam::operator>>(Istream& is, Matrix<T>& M)
+template<class Form, class Type>
+Foam::Istream& Foam::operator>>(Istream& is, Matrix<Form, Type>& M)
 {
     // Anull matrix
     M.clear();
 
-    is.fatalCheck("operator>>(Istream&, Matrix<T>&)");
+    is.fatalCheck("operator>>(Istream&, Matrix<Form, Type>&)");
 
     token firstToken(is);
 
-    is.fatalCheck("operator>>(Istream&, Matrix<T>&) : reading first token");
+    is.fatalCheck("operator>>(Istream&, Matrix<Form, Type>&) : reading first token");
 
     if (firstToken.isLabel())
     {
@@ -63,7 +63,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<T>& M)
         label nm = M.n_*M.m_;
 
         // Read list contents depending on data format
-        if (is.format() == IOstream::ASCII || !contiguous<T>())
+        if (is.format() == IOstream::ASCII || !contiguous<Type>())
         {
             // Read beginning of contents
             char listDelimiter = is.readBeginList("Matrix");
@@ -71,7 +71,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<T>& M)
             if (nm)
             {
                 M.allocate();
-                T* v = M.v_[0];
+                Type* v = M.v_[0];
 
                 if (listDelimiter == token::BEGIN_LIST)
                 {
@@ -88,7 +88,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<T>& M)
 
                             is.fatalCheck
                             (
-                                "operator>>(Istream&, Matrix<T>&) : "
+                                "operator>>(Istream&, Matrix<Form, Type>&) : "
                                 "reading entry"
                             );
                         }
@@ -98,12 +98,12 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<T>& M)
                 }
                 else
                 {
-                    T element;
+                    Type element;
                     is >> element;
 
                     is.fatalCheck
                     (
-                        "operator>>(Istream&, Matrix<T>&) : "
+                        "operator>>(Istream&, Matrix<Form, Type>&) : "
                         "reading the single entry"
                     );
 
@@ -122,13 +122,13 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<T>& M)
             if (nm)
             {
                 M.allocate();
-                T* v = M.v_[0];
+                Type* v = M.v_[0];
 
-                is.read(reinterpret_cast<char*>(v), nm*sizeof(T));
+                is.read(reinterpret_cast<char*>(v), nm*sizeof(Type));
 
                 is.fatalCheck
                 (
-                    "operator>>(Istream&, Matrix<T>&) : "
+                    "operator>>(Istream&, Matrix<Form, Type>&) : "
                     "reading the binary block"
                 );
             }
@@ -136,7 +136,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<T>& M)
     }
     else
     {
-        FatalIOErrorIn("operator>>(Istream&, Matrix<T>&)", is)
+        FatalIOErrorIn("operator>>(Istream&, Matrix<Form, Type>&)", is)
             << "incorrect first token, expected <int>, found "
             << firstToken.info()
             << exit(FatalIOError);
@@ -146,23 +146,23 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<T>& M)
 }
 
 
-template<class T>
-Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<T>& M)
+template<class Form, class Type>
+Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<Form, Type>& M)
 {
     label nm = M.n_*M.m_;
 
     os  << M.n() << token::SPACE << M.m();
 
     // Write list contents depending on data format
-    if (os.format() == IOstream::ASCII || !contiguous<T>())
+    if (os.format() == IOstream::ASCII || !contiguous<Type>())
     {
         if (nm)
         {
             bool uniform = false;
 
-            const T* v = M.v_[0];
+            const Type* v = M.v_[0];
 
-            if (nm > 1 && contiguous<T>())
+            if (nm > 1 && contiguous<Type>())
             {
                 uniform = true;
 
@@ -187,7 +187,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<T>& M)
                 // Write end of contents delimiter
                 os << token::END_BLOCK;
             }
-            else if (nm < 10 && contiguous<T>())
+            else if (nm < 10 && contiguous<Type>())
             {
                 // Write size of list and start contents delimiter
                 os  << token::BEGIN_LIST;
@@ -246,7 +246,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<T>& M)
     {
         if (nm)
         {
-            os.write(reinterpret_cast<const char*>(M.v_[0]), nm*sizeof(T));
+            os.write(reinterpret_cast<const char*>(M.v_[0]), nm*sizeof(Type));
         }
     }
 

src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H

@@ -0,0 +1,92 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::RectangularMatrix
+
+Description
+    A templated 2D rectangular matrix of objects of \<T\>, where the n x n matrix
+    dimension is known and used for subscript bounds checking, etc.
+
+SourceFiles
+    RectangularMatrixI.H
+    RectangularMatrix.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef RectangularMatrix_H
+#define RectangularMatrix_H
+
+#include "Matrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class Matrix Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class RectangularMatrix
+:
+    public Matrix<RectangularMatrix<Type>, Type>
+{
+
+public:
+
+    // Constructors
+
+        //- Null constructor.
+        inline RectangularMatrix();
+
+        //- Construct given number of rows and columns,
+        inline RectangularMatrix(const label m, const label n);
+
+        //- Construct with given number of rows and columns
+        //  and value for all elements.
+        inline RectangularMatrix(const label m, const label n, const Type&);
+
+        //- Construct from Istream.
+        inline RectangularMatrix(Istream&);
+
+        //- Clone
+        inline autoPtr<RectangularMatrix<Type> > clone() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#   include "RectangularMatrixI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H

@@ -0,0 +1,70 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+inline Foam::RectangularMatrix<Type>::RectangularMatrix()
+:
+    Matrix<RectangularMatrix<Type>, Type>()
+{}
+
+template<class Type>
+inline Foam::RectangularMatrix<Type>::RectangularMatrix
+(
+    const label m,
+    const label n
+)
+:
+    Matrix<RectangularMatrix<Type>, Type>(m, n)
+{}
+
+template<class Type>
+inline Foam::RectangularMatrix<Type>::RectangularMatrix
+(
+    const label m,
+    const label n,
+    const Type& t
+)
+:
+    Matrix<RectangularMatrix<Type>, Type>(m, n, t)
+{}
+
+template<class Type>
+inline Foam::RectangularMatrix<Type>::RectangularMatrix(Istream& is)
+:
+    Matrix<RectangularMatrix<Type>, Type>(is)
+{}
+
+template<class Type>
+inline Foam::autoPtr<Foam::RectangularMatrix<Type> >
+Foam::RectangularMatrix<Type>::clone() const
+{
+    return autoPtr<RectangularMatrix<Type> >(new RectangularMatrix<Type>(*this));
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H

@@ -23,22 +23,22 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::scalarMatrix
+    Foam::SquareMatrix
 
 Description
-    Foam::scalarMatrix
+    A templated 2D square matrix of objects of \<T\>, where the n x n matrix
+    dimension is known and used for subscript bounds checking, etc.
 
 SourceFiles
-    scalarMatrix.C
+    SquareMatrixI.H
+    SquareMatrix.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef scalarMatrix_H
-#define scalarMatrix_H
+#ifndef SquareMatrix_H
+#define SquareMatrix_H
 
 #include "Matrix.H"
-#include "scalarField.H"
-#include "labelList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -46,65 +46,39 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class scalarMatrix Declaration
+                           Class Matrix Declaration
 \*---------------------------------------------------------------------------*/
 
-class scalarMatrix
+template<class Type>
+class SquareMatrix
 :
-    public Matrix<scalar>
+    public Matrix<SquareMatrix<Type>, Type>
 {
 
 public:
 
     // Constructors
 
-        //- Construct null
-        scalarMatrix();
+        //- Null constructor.
+        inline SquareMatrix();
 
-        //- Construct given size
-        scalarMatrix(const label);
+        //- Construct given number of rows/columns.
+        inline SquareMatrix(const label n);
 
-        //- Construct from Matrix<scalar>
-        scalarMatrix(const Matrix<scalar>&);
+        //- Construct given number of rows and columns,
+        //  It checks that m == n.
+        inline SquareMatrix(const label m, const label n);
 
-        //- Construct from Istream
-        scalarMatrix(Istream&);
+        //- Construct with given number of rows and rows
+        //  and value for all elements.
+        //  It checks that m == n.
+        inline SquareMatrix(const label m, const label n, const Type&);
 
+        //- Construct from Istream.
+        inline SquareMatrix(Istream&);
 
-    // Member Functions
-
-        //- Solve the matrix using Gaussian elimination with pivoting,
-        //  returning the solution in the source
-        template<class T>
-        static void solve(Matrix<scalar>& matrix, Field<T>& source);
-
-        //- Solve the matrix using Gaussian elimination with pivoting
-        //  and return the solution
-        template<class T>
-        void solve(Field<T>& psi, const Field<T>& source) const;
-
-
-        //- LU decompose the matrix with pivoting
-        static void LUDecompose
-        (
-            Matrix<scalar>& matrix,
-            labelList& pivotIndices
-        );
-
-        //- LU back-substitution with given source, returning the solution
-        //  in the source
-        template<class T>
-        static void LUBacksubstitute
-        (
-            const Matrix<scalar>& luMmatrix,
-            const labelList& pivotIndices,
-            Field<T>& source
-        );
-
-        //- Solve the matrix using LU decomposition with pivoting
-        //  returning the LU form of the matrix and the solution in the source
-        template<class T>
-        static void LUsolve(Matrix<scalar>& matrix, Field<T>& source);
+        //- Clone
+        inline autoPtr<SquareMatrix<Type> > clone() const;
 };
 
 
@@ -114,9 +88,7 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#ifdef NoRepository
-#   include "scalarMatrixTemplates.C"
-#endif
+#   include "SquareMatrixI.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H

@@ -0,0 +1,89 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+inline Foam::SquareMatrix<Type>::SquareMatrix()
+:
+    Matrix<SquareMatrix<Type>, Type>()
+{}
+
+template<class Type>
+inline Foam::SquareMatrix<Type>::SquareMatrix(const label n)
+:
+    Matrix<SquareMatrix<Type>, Type>(n, n)
+{}
+
+template<class Type>
+inline Foam::SquareMatrix<Type>::SquareMatrix(const label m, const label n)
+:
+    Matrix<SquareMatrix<Type>, Type>(m, n)
+{
+    if (m != n)
+    {
+        FatalErrorIn
+        (
+            "SquareMatrix<Type>::SquareMatrix(const label m, const label n)"
+        ) << "m != n for constructing a square matrix" << exit(FatalError);
+    }
+}
+
+template<class Type>
+inline Foam::SquareMatrix<Type>::SquareMatrix
+(
+    const label m,
+    const label n,
+    const Type& t
+)
+:
+    Matrix<SquareMatrix<Type>, Type>(m, n, t)
+{
+    if (m != n)
+    {
+        FatalErrorIn
+        (
+            "SquareMatrix<Type>::SquareMatrix"
+            "(const label m, const label n, const Type&)"
+        ) << "m != n for constructing a square matrix" << exit(FatalError);
+    }
+}
+
+template<class Type>
+inline Foam::SquareMatrix<Type>::SquareMatrix(Istream& is)
+:
+    Matrix<SquareMatrix<Type>, Type>(is)
+{}
+
+template<class Type>
+inline Foam::autoPtr<Foam::SquareMatrix<Type> >
+Foam::SquareMatrix<Type>::clone() const
+{
+    return autoPtr<SquareMatrix<Type> >(new SquareMatrix<Type>(*this));
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C

@@ -0,0 +1,403 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2005 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "SVD.H"
+#include "scalarList.H"
+#include "scalarMatrices.H"
+#include "ListOps.H"
+
+// * * * * * * * * * * * * * * * * Constructor  * * * * * * * * * * * * * * //
+
+Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
+:
+    U_(A),
+    V_(A.m(), A.m()),
+    S_(A.m()),
+    VSinvUt_(A.m(), A.n()),
+    nZeros_(0)
+{
+    // SVDcomp to find U_, V_ and S_ - the singular values
+
+    const label Um = U_.m();
+    const label Un = U_.n();
+
+    scalarList rv1(Um);
+    scalar g = 0;
+    scalar scale = 0;
+    scalar s = 0;
+    scalar anorm = 0;
+    label l = 0;
+
+    for (label i = 0; i < Um; i++)
+    {
+        l = i+2;
+        rv1[i] = scale*g;
+        g = s = scale = 0;
+
+        if (i < Un)
+        {
+            for (label k = i; k < Un; k++)
+            {
+                scale += mag(U_[k][i]);
+            }
+
+            if (scale != 0)
+            {
+                for (label k = i; k < Un; k++)
+                {
+                    U_[k][i] /= scale;
+                    s += U_[k][i]*U_[k][i];
+                }
+
+                scalar f = U_[i][i];
+                g = -sign(Foam::sqrt(s), f);
+                scalar h = f*g - s;
+                U_[i][i] = f - g;
+
+                for (label j = l-1; j < Um; j++)
+                {
+                    s = 0;
+                    for (label k = i; k < Un; k++)
+                    {
+                        s += U_[k][i]*U_[k][j];
+                    }
+
+                    f = s/h;
+                    for (label k = i; k < A.n(); k++)
+                    {
+                        U_[k][j] += f*U_[k][i];
+                    }
+                }
+
+                for (label k = i; k < Un; k++)
+                {
+                    U_[k][i] *= scale;
+                }
+            }
+        }
+
+        S_[i] = scale*g;
+
+        g = s = scale = 0;
+
+        if (i+1 <= Un && i != Um)
+        {
+            for (label k = l-1; k < Um; k++)
+            {
+                scale += mag(U_[i][k]);
+            }
+
+            if (scale != 0)
+            {
+                for (label k=l-1; k < Um; k++)
+                {
+                    U_[i][k] /= scale;
+                    s += U_[i][k]*U_[i][k];
+                }
+
+                scalar f = U_[i][l-1];
+                g = -sign(Foam::sqrt(s),f);
+                scalar h = f*g - s;
+                U_[i][l-1] = f - g;
+
+                for (label k = l-1; k < Um; k++)
+                {
+                    rv1[k] = U_[i][k]/h;
+                }
+
+                for (label j = l-1; j < Un; j++)
+                {
+                    s = 0;
+                    for (label k = l-1; k < Um; k++)
+                    {
+                        s += U_[j][k]*U_[i][k];
+                    }
+
+                    for (label k = l-1; k < Um; k++)
+                    {
+                        U_[j][k] += s*rv1[k];
+                    }
+                }
+                for (label k = l-1; k < Um; k++)
+                {
+                    U_[i][k] *= scale;
+                }
+            }
+        }
+
+        anorm = max(anorm, mag(S_[i]) + mag(rv1[i]));
+    }
+
+    for (label i = Um-1; i >= 0; i--)
+    {
+        if (i < Um-1)
+        {
+            if (g != 0)
+            {
+                for (label j = l; j < Um; j++)
+                {
+                    V_[j][i] = (U_[i][j]/U_[i][l])/g;
+                }
+
+                for (label j=l; j < Um; j++)
+                {
+                    s = 0;
+                    for (label k = l; k < Um; k++)
+                    {
+                        s += U_[i][k]*V_[k][j];
+                    }
+
+                    for (label k = l; k < Um; k++)
+                    {
+                        V_[k][j] += s*V_[k][i];
+                    }
+                }
+            }
+
+            for (label j = l; j < Um;j++)
+            {
+                V_[i][j] = V_[j][i] = 0.0;
+            }
+        }
+
+        V_[i][i] = 1;
+        g = rv1[i];
+        l = i;
+    }
+
+    for (label i = min(Um, Un) - 1; i >= 0; i--)
+    {
+        l = i+1;
+        g = S_[i];
+
+        for (label j = l; j < Um; j++)
+        {
+            U_[i][j] = 0.0;
+        }
+
+        if (g != 0)
+        {
+            g = 1.0/g;
+
+            for (label j = l; j < Um; j++)
+            {
+                s = 0;
+                for (label k = l; k < Un; k++)
+                {
+                    s += U_[k][i]*U_[k][j];
+                }
+
+                scalar f = (s/U_[i][i])*g;
+
+                for (label k = i; k < Un; k++)
+                {
+                    U_[k][j] += f*U_[k][i];
+                }
+            }
+
+            for (label j = i; j < Un; j++)
+            {
+                U_[j][i] *= g;
+            }
+        }
+        else
+        {
+            for (label j = i; j < Un; j++)
+            {
+                U_[j][i] = 0.0;
+            }
+        }
+
+        ++U_[i][i];
+    }
+
+    for (label k = Um-1; k >= 0; k--)
+    {
+        for (label its = 0; its < 35; its++)
+        {
+            bool flag = true;
+
+            label nm;
+            for (l = k; l >= 0; l--)
+            {
+                nm = l-1;
+                if (mag(rv1[l]) + anorm == anorm)
+                {
+                    flag = false;
+                    break;
+                }
+                if (mag(S_[nm]) + anorm == anorm) break;
+            }
+
+            if (flag)
+            {
+                scalar c = 0.0;
+                s = 1.0;
+                for (label i = l-1; i < k+1; i++)
+                {
+                    scalar f = s*rv1[i];
+                    rv1[i] = c*rv1[i];
+
+                    if (mag(f) + anorm == anorm) break;
+
+                    g = S_[i];
+                    scalar h = sqrtSumSqr(f, g);
+                    S_[i] = h;
+                    h = 1.0/h;
+                    c = g*h;
+                    s = -f*h;
+
+                    for (label j = 0; j < Un; j++)
+                    {
+                        scalar y = U_[j][nm];
+                        scalar z = U_[j][i];
+                        U_[j][nm] = y*c + z*s;
+                        U_[j][i] = z*c - y*s;
+                    }
+                }
+            }
+
+            scalar z = S_[k];
+
+            if (l == k)
+            {
+                if (z < 0.0)
+                {
+                    S_[k] = -z;
+
+                    for (label j = 0; j < Um; j++) V_[j][k] = -V_[j][k];
+                }
+                break;
+            }
+            if (its == 34)
+            {
+                WarningIn
+                (
+                    "SVD::SVD"
+                    "(scalarRectangularMatrix& A, const scalar minCondition)"
+                )   << "no convergence in 35 SVD iterations"
+                    << endl;
+            }
+
+            scalar x = S_[l];
+            nm = k-1;
+            scalar y = S_[nm];
+            g = rv1[nm];
+            scalar h = rv1[k];
+            scalar f = ((y - z)*(y + z) + (g - h)*(g + h))/(2.0*h*y);
+            g = sqrtSumSqr(f, 1.0);
+            f = ((x - z)*(x + z) + h*((y/(f + sign(g, f))) - h))/x;
+            scalar c = 1.0;
+            s = 1.0;
+
+            for (label j = l; j <= nm; j++)
+            {
+                label i = j + 1;
+                g = rv1[i];
+                y = S_[i];
+                h = s*g;
+                g = c*g;
+                scalar z = sqrtSumSqr(f, h);
+                rv1[j] = z;
+                c = f/z;
+                s = h/z;
+                f = x*c + g*s;
+                g = g*c - x*s;
+                h = y*s;
+                y *= c;
+
+                for (label jj = 0; jj < Um; jj++)
+                {
+                    x = V_[jj][j];
+                    z = V_[jj][i];
+                    V_[jj][j] = x*c + z*s;
+                    V_[jj][i] = z*c - x*s;
+                }
+
+                z = sqrtSumSqr(f, h);
+                S_[j] = z;
+                if (z)
+                {
+                    z = 1.0/z;
+                    c = f*z;
+                    s = h*z;
+                }
+                f = c*g + s*y;
+                x = c*y - s*g;
+
+                for (label jj=0; jj < Un; jj++)
+                {
+                    y = U_[jj][j];
+                    z = U_[jj][i];
+                    U_[jj][j] = y*c + z*s;
+                    U_[jj][i] = z*c - y*s;
+                }
+            }
+            rv1[l] = 0.0;
+            rv1[k] = f;
+            S_[k] = x;
+        }
+    }
+
+    // zero singular values that are less than minCondition*maxS
+    const scalar minS = minCondition*S_[findMax(S_)];
+    for (label i = 0; i < S_.size(); i++)
+    {
+        if (S_[i] <= minS)
+        {
+            //Info << "Removing " << S_[i] << " < " << minS << endl;
+            S_[i] = 0;
+            nZeros_++;
+        }
+    }
+
+    // now multiply out to find the pseudo inverse of A, VSinvUt_
+    multiply(VSinvUt_, V_, inv(S_), U_.T());
+
+    // test SVD
+    /*scalarRectangularMatrix SVDA(A.n(), A.m());
+    multiply(SVDA, U_, S_, transpose(V_));
+    scalar maxDiff = 0;
+    scalar diff = 0;
+    for(label i = 0; i < A.n(); i++)
+    {
+        for(label j = 0; j < A.m(); j++)
+        {
+            diff = mag(A[i][j] - SVDA[i][j]);
+            if (diff > maxDiff) maxDiff = diff;
+        }
+    }
+    Info << "Maximum discrepancy between A and svd(A) = " << maxDiff << endl;
+
+    if (maxDiff > 4)
+    {
+        Info << "singular values " << S_ << endl;
+    }
+    */
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.H

@@ -0,0 +1,128 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2005 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::SVD
+
+Description
+    Singular value decomposition of a rectangular matrix.
+
+SourceFiles
+    SVDI.H
+    SVD.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef SVD_H
+#define SVD_H
+
+#include "scalarMatrices.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+/*---------------------------------------------------------------------------*\
+                      Class SVD Declaration
+\*---------------------------------------------------------------------------*/
+
+class SVD
+{
+    // Private data
+
+        //- Rectangular matrix with the same dimensions as the input
+        scalarRectangularMatrix U_;
+
+        //- square matrix V
+        scalarRectangularMatrix V_;
+
+        //- The singular values
+        DiagonalMatrix<scalar> S_;
+
+        //- The matrix product V S^(-1) U^T
+        scalarRectangularMatrix VSinvUt_;
+
+        //- The number of zero singular values
+        label nZeros_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        SVD(const SVD&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const SVD&);
+
+        template<class T>
+        inline const T sign(const T& a, const T& b);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from a rectangular Matrix
+        SVD(const scalarRectangularMatrix& A, const scalar minCondition = 0);
+
+
+    // Access functions
+
+        //- Return U
+        inline const scalarRectangularMatrix& U() const;
+
+        //- Return the square matrix V
+        inline const scalarRectangularMatrix& V() const;
+
+        //- Return the singular values
+        inline const scalarDiagonalMatrix& S() const;
+
+        //- Return VSinvUt (the pseudo inverse)
+        inline const scalarRectangularMatrix& VSinvUt() const;
+
+        //- Return the number of zero singular values
+        inline label nZeros() const;
+
+        //- Return the minimum non-zero singular value
+        inline scalar minNonZeroS() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "SVDI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/matrices/scalarMatrices/SVD/SVDI.H

@@ -0,0 +1,75 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * * //
+
+template<class T>
+inline const T Foam::SVD::sign(const T& a, const T& b)
+{
+    return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::scalarRectangularMatrix& Foam::SVD::U() const
+{
+    return U_;
+}
+
+inline const Foam::scalarRectangularMatrix& Foam::SVD::V() const
+{
+    return V_;
+}
+
+inline const Foam::scalarDiagonalMatrix& Foam::SVD::S() const
+{
+    return S_;
+}
+
+inline const Foam::scalarRectangularMatrix& Foam::SVD::VSinvUt() const
+{
+    return VSinvUt_;
+}
+
+inline Foam::label Foam::SVD::nZeros() const
+{
+    return nZeros_;
+}
+
+inline Foam::scalar Foam::SVD::minNonZeroS() const
+{
+    scalar minS = S_[0];
+    for(label i = 1; i < S_.size(); i++)
+    {
+        scalar s = S_[i];
+        if (s > VSMALL && s < minS) minS = s;
+    }
+    return minS;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C

@@ -0,0 +1,293 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "scalarMatrices.H"
+#include "SVD.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::LUDecompose
+(
+    scalarSquareMatrix& matrix,
+    labelList& pivotIndices
+)
+{
+    label n = matrix.n();
+    scalar vv[n];
+
+    for (register label i=0; i<n; i++)
+    {
+        scalar largestCoeff = 0.0;
+        scalar temp;
+        const scalar* __restrict__ matrixi = matrix[i];
+
+        for (register label j=0; j<n; j++)
+        {
+            if ((temp = mag(matrixi[j])) > largestCoeff)
+            {
+                largestCoeff = temp;
+            }
+        }
+
+        if (largestCoeff == 0.0)
+        {
+            FatalErrorIn
+            (
+                "LUdecompose"
+                "(scalarSquareMatrix& matrix, labelList& rowIndices)"
+            )   << "Singular matrix" << exit(FatalError);
+        }
+
+        vv[i] = 1.0/largestCoeff;
+    }
+
+    for (register label j=0; j<n; j++)
+    {
+        scalar* __restrict__ matrixj = matrix[j];
+
+        for (register label i=0; i<j; i++)
+        {
+            scalar* __restrict__ matrixi = matrix[i];
+
+            scalar sum = matrixi[j];
+            for (register label k=0; k<i; k++)
+            {
+                sum -= matrixi[k]*matrix[k][j];
+            }
+            matrixi[j] = sum;
+        }
+
+        label iMax = 0;
+
+        scalar largestCoeff = 0.0;
+        for (register label i=j; i<n; i++)
+        {
+            scalar* __restrict__ matrixi = matrix[i];
+            scalar sum = matrixi[j];
+
+            for (register label k=0; k<j; k++)
+            {
+                sum -= matrixi[k]*matrix[k][j];
+            }
+
+            matrixi[j] = sum;
+
+            scalar temp;
+            if ((temp = vv[i]*mag(sum)) >= largestCoeff)
+            {
+                largestCoeff = temp;
+                iMax = i;
+            }
+        }
+
+        pivotIndices[j] = iMax;
+
+        if (j != iMax)
+        {
+            scalar* __restrict__ matrixiMax = matrix[iMax];
+
+            for (register label k=0; k<n; k++)
+            {
+                Swap(matrixj[k], matrixiMax[k]);
+            }
+
+            vv[iMax] = vv[j];
+        }
+
+        if (matrixj[j] == 0.0)
+        {
+            matrixj[j] = SMALL;
+        }
+
+        if (j != n-1)
+        {
+            scalar rDiag = 1.0/matrixj[j];
+
+            for (register label i=j+1; i<n; i++)
+            {
+                matrix[i][j] *= rDiag;
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
+
+void Foam::multiply
+(
+    scalarRectangularMatrix& ans,         // value changed in return
+    const scalarRectangularMatrix& A,
+    const scalarRectangularMatrix& B
+)
+{
+    if (A.m() != B.n())
+    {
+        FatalErrorIn
+        (
+            "multiply("
+            "scalarRectangularMatrix& answer "
+            "const scalarRectangularMatrix& A, "
+            "const scalarRectangularMatrix& B)"
+        )   << "A and B must have identical inner dimensions but A.m = "
+            << A.m() << " and B.n = " << B.n()
+            << abort(FatalError);
+    }
+
+    ans = scalarRectangularMatrix(A.n(), B.m(), scalar(0));
+
+    for(register label i = 0; i < A.n(); i++)
+    {
+        for(register label j = 0; j < B.m(); j++)
+        {
+            for(register label l = 0; l < B.n(); l++)
+            {
+                ans[i][j] += A[i][l]*B[l][j];
+            }
+        }
+    }
+}
+
+
+void Foam::multiply
+(
+    scalarRectangularMatrix& ans,         // value changed in return
+    const scalarRectangularMatrix& A,
+    const scalarRectangularMatrix& B,
+    const scalarRectangularMatrix& C
+)
+{
+    if (A.m() != B.n())
+    {
+        FatalErrorIn
+        (
+            "multiply("
+            "const scalarRectangularMatrix& A, "
+            "const scalarRectangularMatrix& B, "
+            "const scalarRectangularMatrix& C, "
+            "scalarRectangularMatrix& answer)"
+        )   << "A and B must have identical inner dimensions but A.m = "
+            << A.m() << " and B.n = " << B.n()
+            << abort(FatalError);
+    }
+
+    if (B.m() != C.n())
+    {
+        FatalErrorIn
+        (
+            "multiply("
+            "const scalarRectangularMatrix& A, "
+            "const scalarRectangularMatrix& B, "
+            "const scalarRectangularMatrix& C, "
+            "scalarRectangularMatrix& answer)"
+        )   << "B and C must have identical inner dimensions but B.m = "
+            << B.m() << " and C.n = " << C.n()
+            << abort(FatalError);
+    }
+
+    ans = scalarRectangularMatrix(A.n(), C.m(), scalar(0));
+
+    for(register label i = 0; i < A.n(); i++)
+    {
+        for(register label g = 0; g < C.m(); g++)
+        {
+            for(register label l = 0; l < C.n(); l++)
+            {
+                scalar ab = 0;
+                for(register label j = 0; j < A.m(); j++)
+                {
+                    ab += A[i][j]*B[j][l];
+                }
+                ans[i][g] += C[l][g] * ab;
+            }
+        }
+    }
+}
+
+
+void Foam::multiply
+(
+    scalarRectangularMatrix& ans,         // value changed in return
+    const scalarRectangularMatrix& A,
+    const DiagonalMatrix<scalar>& B,
+    const scalarRectangularMatrix& C
+)
+{
+    if (A.m() != B.size())
+    {
+        FatalErrorIn
+        (
+            "multiply("
+            "const scalarRectangularMatrix& A, "
+            "const DiagonalMatrix<scalar>& B, "
+            "const scalarRectangularMatrix& C, "
+            "scalarRectangularMatrix& answer)"
+        )   << "A and B must have identical inner dimensions but A.m = "
+            << A.m() << " and B.n = " << B.size()
+            << abort(FatalError);
+    }
+
+    if (B.size() != C.n())
+    {
+        FatalErrorIn
+        (
+            "multiply("
+            "const scalarRectangularMatrix& A, "
+            "const DiagonalMatrix<scalar>& B, "
+            "const scalarRectangularMatrix& C, "
+            "scalarRectangularMatrix& answer)"
+        )   << "B and C must have identical inner dimensions but B.m = "
+            << B.size() << " and C.n = " << C.n()
+            << abort(FatalError);
+    }
+
+    ans = scalarRectangularMatrix(A.n(), C.m(), scalar(0));
+
+    for(register label i = 0; i < A.n(); i++)
+    {
+        for(register label g = 0; g < C.m(); g++)
+        {
+            for(register label l = 0; l < C.n(); l++)
+            {
+                ans[i][g] += C[l][g] * A[i][l]*B[l];
+            }
+        }
+    }
+}
+
+
+Foam::RectangularMatrix<Foam::scalar> Foam::SVDinv
+(
+    const scalarRectangularMatrix& A,
+    scalar minCondition
+)
+{
+    SVD svd(A, minCondition);
+    return svd.VSinvUt();
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H

@@ -0,0 +1,137 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    scalarMatrices
+
+Description
+    Scalar matrices
+
+SourceFiles
+    scalarMatrices.C
+    scalarMatricesTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef scalarMatrices_H
+#define scalarMatrices_H
+
+#include "RectangularMatrix.H"
+#include "SquareMatrix.H"
+#include "DiagonalMatrix.H"
+#include "scalarField.H"
+#include "labelList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+typedef RectangularMatrix<scalar> scalarRectangularMatrix;
+typedef SquareMatrix<scalar> scalarSquareMatrix;
+typedef DiagonalMatrix<scalar> scalarDiagonalMatrix;
+
+//- Solve the matrix using Gaussian elimination with pivoting,
+//  returning the solution in the source
+template<class Type>
+void solve(scalarSquareMatrix& matrix, Field<Type>& source);
+
+//- Solve the matrix using Gaussian elimination with pivoting
+//  and return the solution
+template<class Type>
+void solve
+(
+    Field<Type>& psi,
+    const scalarSquareMatrix& matrix,
+    const Field<Type>& source
+);
+
+//- LU decompose the matrix with pivoting
+void LUDecompose
+(
+    scalarSquareMatrix& matrix,
+    labelList& pivotIndices
+);
+
+//- LU back-substitution with given source, returning the solution
+//  in the source
+template<class Type>
+void LUBacksubstitute
+(
+    const scalarSquareMatrix& luMmatrix,
+    const labelList& pivotIndices,
+    Field<Type>& source
+);
+
+//- Solve the matrix using LU decomposition with pivoting
+//  returning the LU form of the matrix and the solution in the source
+template<class Type>
+void LUsolve(scalarSquareMatrix& matrix, Field<Type>& source);
+
+void multiply
+(
+    scalarRectangularMatrix& answer,         // value changed in return
+    const scalarRectangularMatrix& A,
+    const scalarRectangularMatrix& B
+);
+
+void multiply
+(
+    scalarRectangularMatrix& answer,         // value changed in return
+    const scalarRectangularMatrix& A,
+    const scalarRectangularMatrix& B,
+    const scalarRectangularMatrix& C
+);
+
+void multiply
+(
+    scalarRectangularMatrix& answer,         // value changed in return
+    const scalarRectangularMatrix& A,
+    const DiagonalMatrix<scalar>& B,
+    const scalarRectangularMatrix& C
+);
+
+//- Return the inverse of matrix A using SVD
+scalarRectangularMatrix SVDinv
+(
+    const scalarRectangularMatrix& A,
+    scalar minCondition = 0
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "scalarMatricesTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C

@@ -24,16 +24,16 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "scalarMatrix.H"
+#include "scalarMatrices.H"
 #include "Swap.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class T>
-void Foam::scalarMatrix::solve
+template<class Type>
+void Foam::solve
 (
-    Matrix<scalar>& tmpMatrix,
-    Field<T>& sourceSol
+    scalarSquareMatrix& tmpMatrix,
+    Field<Type>& sourceSol
 )
 {
     label n = tmpMatrix.n();
@@ -68,7 +68,7 @@ void Foam::scalarMatrix::solve
         // Check that the system of equations isn't singular
         if (mag(tmpMatrix[i][i]) < 1e-20)
         {
-            FatalErrorIn("scalarMatrix::solve()")
+            FatalErrorIn("solve(scalarSquareMatrix&, Field<Type>& sourceSol)")
                 << "Singular Matrix"
                 << exit(FatalError);
         }
@@ -89,7 +89,7 @@ void Foam::scalarMatrix::solve
     // Back-substitution
     for (register label j=n-1; j>=0; j--)
     {
-        T ntempvec = pTraits<T>::zero;
+        Type ntempvec = pTraits<Type>::zero;
 
         for (register label k=j+1; k<n; k++)
         {
@@ -101,21 +101,26 @@ void Foam::scalarMatrix::solve
 }
 
 
-template<class T>
-void Foam::scalarMatrix::solve(Field<T>& psi, const Field<T>& source) const
+template<class Type>
+void Foam::solve
+(
+    Field<Type>& psi,
+    const scalarSquareMatrix& matrix,
+    const Field<Type>& source
+)
 {
-    Matrix<scalar> tmpMatrix = *this;
+    scalarSquareMatrix tmpMatrix = matrix;
     psi = source;
     solve(tmpMatrix, psi);
 }
 
 
-template<class T>
-void Foam::scalarMatrix::LUBacksubstitute
+template<class Type>
+void Foam::LUBacksubstitute
 (
-    const Matrix<scalar>& luMatrix,
+    const scalarSquareMatrix& luMatrix,
     const labelList& pivotIndices,
-    Field<T>& sourceSol
+    Field<Type>& sourceSol
 )
 {
     label n = luMatrix.n();
@@ -125,7 +130,7 @@ void Foam::scalarMatrix::LUBacksubstitute
     for (register label i=0; i<n; i++)
     {
         label ip = pivotIndices[i];
-        T sum = sourceSol[ip];
+        Type sum = sourceSol[ip];
         sourceSol[ip] = sourceSol[i];
         const scalar* __restrict__ luMatrixi = luMatrix[i];
 
@@ -136,7 +141,7 @@ void Foam::scalarMatrix::LUBacksubstitute
                 sum -= luMatrixi[j]*sourceSol[j];
             }
         }
-        else if (sum != pTraits<T>::zero)
+        else if (sum != pTraits<Type>::zero)
         {
             ii = i+1;
         }
@@ -146,11 +151,11 @@ void Foam::scalarMatrix::LUBacksubstitute
 
     for (register label i=n-1; i>=0; i--)
     {
-        T sum = sourceSol[i];
+        Type sum = sourceSol[i];
         const scalar* __restrict__ luMatrixi = luMatrix[i];
 
         for (register label j=i+1; j<n; j++)
-        { 
+        {
             sum -= luMatrixi[j]*sourceSol[j];
         }
 
@@ -159,11 +164,11 @@ void Foam::scalarMatrix::LUBacksubstitute
 }
 
 
-template<class T>
-void Foam::scalarMatrix::LUsolve
+template<class Type>
+void Foam::LUsolve
 (
-    Matrix<scalar>& matrix,
-    Field<T>& sourceSol
+    scalarSquareMatrix& matrix,
+    Field<Type>& sourceSol
 )
 {
     labelList pivotIndices(matrix.n());

src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.C

@@ -1,161 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "scalarMatrix.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::scalarMatrix::scalarMatrix()
-{}
-
-
-Foam::scalarMatrix::scalarMatrix(const label mSize)
-:
-    Matrix<scalar>(mSize, mSize, 0.0)
-{}
-
-
-Foam::scalarMatrix::scalarMatrix(const Matrix<scalar>& matrix)
-:
-    Matrix<scalar>(matrix)
-{}
-
-
-Foam::scalarMatrix::scalarMatrix(Istream& is)
-:
-    Matrix<scalar>(is)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::scalarMatrix::LUDecompose
-(
-    Matrix<scalar>& matrix,
-    labelList& pivotIndices
-)
-{
-    label n = matrix.n();
-    scalar vv[n];
-
-    for (register label i=0; i<n; i++)
-    {
-        scalar largestCoeff = 0.0;
-        scalar temp;
-        const scalar* __restrict__ matrixi = matrix[i];
-
-        for (register label j=0; j<n; j++)
-        {
-            if ((temp = mag(matrixi[j])) > largestCoeff)
-            {
-                largestCoeff = temp;
-            }
-        }
-
-        if (largestCoeff == 0.0)
-        {
-            FatalErrorIn
-            (
-                "scalarMatrix::LUdecompose"
-                "(Matrix<scalar>& matrix, labelList& rowIndices)"
-            )   << "Singular matrix" << exit(FatalError);
-        }
-
-        vv[i] = 1.0/largestCoeff;
-    }
-
-    for (register label j=0; j<n; j++)
-    {
-        scalar* __restrict__ matrixj = matrix[j];
-
-        for (register label i=0; i<j; i++)
-        {
-            scalar* __restrict__ matrixi = matrix[i];
-
-            scalar sum = matrixi[j];
-            for (register label k=0; k<i; k++)
-            {
-                sum -= matrixi[k]*matrix[k][j];
-            }
-            matrixi[j] = sum;
-        }
-
-        label iMax = 0;
-
-        scalar largestCoeff = 0.0;
-        for (register label i=j; i<n; i++)
-        {
-            scalar* __restrict__ matrixi = matrix[i];
-            scalar sum = matrixi[j];
-
-            for (register label k=0; k<j; k++)
-            {
-                sum -= matrixi[k]*matrix[k][j];
-            }
-
-            matrixi[j] = sum;
-
-            scalar temp;
-            if ((temp = vv[i]*mag(sum)) >= largestCoeff)
-            {
-                largestCoeff = temp;
-                iMax = i;
-            }
-        }
-
-        pivotIndices[j] = iMax;
-
-        if (j != iMax)
-        {
-            scalar* __restrict__ matrixiMax = matrix[iMax];
-
-            for (register label k=0; k<n; k++)
-            {
-                Swap(matrixj[k], matrixiMax[k]);
-            }
-            
-            vv[iMax] = vv[j];
-        }
-
-        if (matrixj[j] == 0.0)
-        {
-            matrixj[j] = SMALL;
-        }
-
-        if (j != n-1)
-        {
-            scalar rDiag = 1.0/matrixj[j];
-
-            for (register label i=j+1; i<n; i++)
-            {
-                matrix[i][j] *= rDiag;
-            }
-        }
-    }
-}
-
-
-// ************************************************************************* //

src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C

@@ -28,55 +28,55 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-template<class T>
-Foam::simpleMatrix<T>::simpleMatrix(const label mSize)
+template<class Type>
+Foam::simpleMatrix<Type>::simpleMatrix(const label mSize)
 :
-    scalarMatrix(mSize),
-    source_(mSize, pTraits<T>::zero)
+    scalarSquareMatrix(mSize),
+    source_(mSize, pTraits<Type>::zero)
 {}
 
 
-template<class T>
-Foam::simpleMatrix<T>::simpleMatrix
+template<class Type>
+Foam::simpleMatrix<Type>::simpleMatrix
 (
-    const scalarMatrix& matrix,
-    const Field<T>& source
+    const scalarSquareMatrix& matrix,
+    const Field<Type>& source
 )
 :
-    scalarMatrix(matrix),
+    scalarSquareMatrix(matrix),
     source_(source)
 {}
 
 
-template<class T>
-Foam::simpleMatrix<T>::simpleMatrix(Istream& is)
+template<class Type>
+Foam::simpleMatrix<Type>::simpleMatrix(Istream& is)
 :
-    scalarMatrix(is),
+    scalarSquareMatrix(is),
     source_(is)
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class T>
-Foam::Field<T> Foam::simpleMatrix<T>::solve() const
+template<class Type>
+Foam::Field<Type> Foam::simpleMatrix<Type>::solve() const
 {
-    scalarMatrix tmpMatrix = *this;
-    Field<T> sourceSol = source_;
+    scalarSquareMatrix tmpMatrix = *this;
+    Field<Type> sourceSol = source_;
 
-    scalarMatrix::solve(tmpMatrix, sourceSol);
+    Foam::solve(tmpMatrix, sourceSol);
 
     return sourceSol;
 }
 
 
-template<class T>
-Foam::Field<T> Foam::simpleMatrix<T>::LUsolve() const
+template<class Type>
+Foam::Field<Type> Foam::simpleMatrix<Type>::LUsolve() const
 {
-    scalarMatrix luMatrix = *this;
-    Field<T> sourceSol = source_;
+    scalarSquareMatrix luMatrix = *this;
+    Field<Type> sourceSol = source_;
 
-    scalarMatrix::LUsolve(luMatrix, sourceSol);
+    Foam::LUsolve(luMatrix, sourceSol);
 
     return sourceSol;
 }
@@ -84,82 +84,82 @@ Foam::Field<T> Foam::simpleMatrix<T>::LUsolve() const
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
-template<class T>
-void Foam::simpleMatrix<T>::operator=(const simpleMatrix<T>& m)
+template<class Type>
+void Foam::simpleMatrix<Type>::operator=(const simpleMatrix<Type>& m)
 {
     if (this == &m)
     {
-        FatalErrorIn("simpleMatrix<T>::operator=(const simpleMatrix<T>&)")
+        FatalErrorIn("simpleMatrix<Type>::operator=(const simpleMatrix<Type>&)")
             << "Attempted assignment to self"
             << abort(FatalError);
     }
 
     if (n() != m.n())
     {
-        FatalErrorIn("simpleMatrix<T>::operator=(const simpleMatrix<T>&)")
+        FatalErrorIn("simpleMatrix<Type>::operator=(const simpleMatrix<Type>&)")
             << "Different size matrices"
             << abort(FatalError);
     }
 
     if (source_.size() != m.source_.size())
     {
-        FatalErrorIn("simpleMatrix<T>::operator=(const simpleMatrix<T>&)")
+        FatalErrorIn("simpleMatrix<Type>::operator=(const simpleMatrix<Type>&)")
             << "Different size source vectors"
             << abort(FatalError);
     }
 
-    scalarMatrix::operator=(m);
+    scalarSquareMatrix::operator=(m);
     source_ = m.source_;
 }
 
 
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
-template<class T>
-Foam::simpleMatrix<T> Foam::operator+
+template<class Type>
+Foam::simpleMatrix<Type> Foam::operator+
 (
-    const simpleMatrix<T>& m1,
-    const simpleMatrix<T>& m2
+    const simpleMatrix<Type>& m1,
+    const simpleMatrix<Type>& m2
 )
 {
-    return simpleMatrix<T>
+    return simpleMatrix<Type>
     (
-        static_cast<const scalarMatrix&>(m1)
-      + static_cast<const scalarMatrix&>(m2),
+        static_cast<const scalarSquareMatrix&>(m1)
+      + static_cast<const scalarSquareMatrix&>(m2),
         m1.source_ + m2.source_
     );
 }
 
 
-template<class T>
-Foam::simpleMatrix<T> Foam::operator-
+template<class Type>
+Foam::simpleMatrix<Type> Foam::operator-
 (
-    const simpleMatrix<T>& m1,
-    const simpleMatrix<T>& m2
+    const simpleMatrix<Type>& m1,
+    const simpleMatrix<Type>& m2
 )
 {
-    return simpleMatrix<T>
+    return simpleMatrix<Type>
     (
-        static_cast<const scalarMatrix&>(m1)
-      - static_cast<const scalarMatrix&>(m2),
+        static_cast<const scalarSquareMatrix&>(m1)
+      - static_cast<const scalarSquareMatrix&>(m2),
         m1.source_ - m2.source_
     );
 }
 
 
-template<class T>
-Foam::simpleMatrix<T> Foam::operator*(const scalar s, const simpleMatrix<T>& m)
+template<class Type>
+Foam::simpleMatrix<Type> Foam::operator*(const scalar s, const simpleMatrix<Type>& m)
 {
-    return simpleMatrix<T>(s*m.matrix_, s*m.source_);
+    return simpleMatrix<Type>(s*m.matrix_, s*m.source_);
 }
 
 
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
-template<class T>
-Foam::Ostream& Foam::operator<<(Ostream& os, const simpleMatrix<T>& m)
+template<class Type>
+Foam::Ostream& Foam::operator<<(Ostream& os, const simpleMatrix<Type>& m)
 {
-    os << static_cast<const scalarMatrix&>(m) << nl << m.source_;
+    os << static_cast<const scalarSquareMatrix&>(m) << nl << m.source_;
     return os;
 }
 

src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H

@@ -36,7 +36,7 @@ SourceFiles
 #ifndef simpleMatrix_H
 #define simpleMatrix_H
 
-#include "scalarMatrix.H"
+#include "scalarMatrices.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -45,35 +45,14 @@ namespace Foam
 
 // Forward declaration of friend functions and operators
 
-template<class T>
+template<class Type>
 class simpleMatrix;
 
-template<class T>
-simpleMatrix<T> operator+
-(
-    const simpleMatrix<T>&,
-    const simpleMatrix<T>&
-);
-
-template<class T>
-simpleMatrix<T> operator-
-(
-    const simpleMatrix<T>&,
-    const simpleMatrix<T>&
-);
-
-template<class T>
-simpleMatrix<T> operator*
-(
-    const scalar,
-    const simpleMatrix<T>&
-);
-
-template<class T>
+template<class Type>
 Ostream& operator<<
 (
     Ostream&,
-    const simpleMatrix<T>&
+    const simpleMatrix<Type>&
 );
 
 
@@ -81,14 +60,14 @@ Ostream& operator<<
                            Class simpleMatrix Declaration
 \*---------------------------------------------------------------------------*/
 
-template<class T>
+template<class Type>
 class simpleMatrix
 :
-    public scalarMatrix
+    public scalarSquareMatrix
 {
     // Private data
 
-        Field<T> source_;
+        Field<Type> source_;
 
 
 public:
@@ -99,25 +78,25 @@ public:
         simpleMatrix(const label);
 
         //- Construct from components
-        simpleMatrix(const scalarMatrix&, const Field<T>&);
+        simpleMatrix(const scalarSquareMatrix&, const Field<Type>&);
 
         //- Construct from Istream
         simpleMatrix(Istream&);
 
         //- Construct as copy
-        simpleMatrix(const simpleMatrix<T>&);
+        simpleMatrix(const simpleMatrix<Type>&);
 
 
     // Member Functions
 
         // Access
 
-            Field<T>& source()
+            Field<Type>& source()
             {
                 return source_;
             }
 
-            const Field<T>& source() const
+            const Field<Type>& source() const
             {
                 return source_;
             }
@@ -125,49 +104,52 @@ public:
 
         //- Solve the matrix using Gaussian elimination with pivoting
         //  and return the solution
-        Field<T> solve() const;
+        Field<Type> solve() const;
 
         //- Solve the matrix using LU decomposition with pivoting
         //  and return the solution
-        Field<T> LUsolve() const;
+        Field<Type> LUsolve() const;
 
 
     // Member Operators
 
-        void operator=(const simpleMatrix<T>&);
-
-
-    // Friend Operators
-
-        friend simpleMatrix<T> operator+ <T>
-        (
-            const simpleMatrix<T>&,
-            const simpleMatrix<T>&
-        );
-
-        friend simpleMatrix<T> operator- <T>
-        (
-            const simpleMatrix<T>&,
-            const simpleMatrix<T>&
-        );
-
-        friend simpleMatrix<T> operator* <T>
-        (
-            const scalar,
-            const simpleMatrix<T>&
-        );
+        void operator=(const simpleMatrix<Type>&);
 
 
     // Ostream Operator
 
-        friend Ostream& operator<< <T>
+        friend Ostream& operator<< <Type>
         (
             Ostream&,
-            const simpleMatrix<T>&
+            const simpleMatrix<Type>&
         );
 };
 
 
+// Global operators
+
+template<class Type>
+simpleMatrix<Type> operator+
+(
+    const simpleMatrix<Type>&,
+    const simpleMatrix<Type>&
+);
+
+template<class Type>
+simpleMatrix<Type> operator-
+(
+    const simpleMatrix<Type>&,
+    const simpleMatrix<Type>&
+);
+
+template<class Type>
+simpleMatrix<Type> operator*
+(
+    const scalar,
+    const simpleMatrix<Type>&
+);
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/finiteVolume/Make/files

@@ -38,6 +38,9 @@ fvMeshMapper = fvMesh/fvMeshMapper
 $(fvMeshMapper)/fvPatchMapper.C
 $(fvMeshMapper)/fvSurfaceMapper.C
 
+extendedStencil = fvMesh/extendedStencil
+$(extendedStencil)/extendedStencil.C
+
 fvPatchFields = fields/fvPatchFields
 $(fvPatchFields)/fvPatchField/fvPatchFields.C
 
@@ -165,6 +168,8 @@ $(schemes)/harmonic/harmonic.C
 $(schemes)/localBlended/localBlended.C
 $(schemes)/localMax/localMax.C
 $(schemes)/localMin/localMin.C
+$(schemes)/quadraticFit/quadraticFit.C
+$(schemes)/quadraticFit/quadraticFitData.C
 
 limitedSchemes = $(surfaceInterpolation)/limitedSchemes
 $(limitedSchemes)/limitedSurfaceInterpolationScheme/limitedSurfaceInterpolationSchemes.C

src/finiteVolume/fvMesh/extendedStencil/extendedStencil.C

@@ -0,0 +1,525 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "extendedStencil.H"
+#include "globalIndex.H"
+#include "syncTools.H"
+#include "SortableList.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Calculates per face a list of global cell/face indices.
+void Foam::extendedStencil::calcFaceStencils
+(
+    const polyMesh& mesh,
+    const globalIndex& globalNumbering
+)
+{
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+    const label nBnd = mesh.nFaces()-mesh.nInternalFaces();
+    const labelList& own = mesh.faceOwner();
+    const labelList& nei = mesh.faceNeighbour();
+
+
+    // Determine neighbouring global cell or boundary face
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    labelList neiGlobal(nBnd);
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+        label faceI = pp.start();
+
+        if (pp.coupled())
+        {
+            // For coupled faces get the cell on the other side
+            forAll(pp, i)
+            {
+                label bFaceI = faceI-mesh.nInternalFaces();
+                neiGlobal[bFaceI] = globalNumbering.toGlobal(own[faceI]);
+                faceI++;
+            }
+        }
+        else if (isA<emptyPolyPatch>(pp))
+        {
+            forAll(pp, i)
+            {
+                label bFaceI = faceI-mesh.nInternalFaces();
+                neiGlobal[bFaceI] = -1;
+                faceI++;
+            }
+        }
+        else
+        {
+            // For noncoupled faces get the boundary face.
+            forAll(pp, i)
+            {
+                label bFaceI = faceI-mesh.nInternalFaces();
+                neiGlobal[bFaceI] =
+                    globalNumbering.toGlobal(mesh.nCells()+bFaceI);
+                faceI++;
+            }
+        }
+    }
+    syncTools::swapBoundaryFaceList(mesh, neiGlobal, false);
+
+
+    // Determine cellCells in global numbering
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    labelListList globalCellCells(mesh.nCells());
+    forAll(globalCellCells, cellI)
+    {
+        const cell& cFaces = mesh.cells()[cellI];
+
+        labelList& cCells = globalCellCells[cellI];
+
+        cCells.setSize(cFaces.size());
+
+        // Collect neighbouring cells/faces
+        label nNbr = 0;
+        forAll(cFaces, i)
+        {
+            label faceI = cFaces[i];
+
+            if (mesh.isInternalFace(faceI))
+            {
+                label nbrCellI = own[faceI];
+                if (nbrCellI == cellI)
+                {
+                    nbrCellI = nei[faceI];
+                }
+                cCells[nNbr++] = globalNumbering.toGlobal(nbrCellI);
+            }
+            else
+            {
+                label nbrCellI = neiGlobal[faceI-mesh.nInternalFaces()];
+                if (nbrCellI != -1)
+                {
+                    cCells[nNbr++] = nbrCellI;
+                }
+            }
+        }
+        cCells.setSize(nNbr);
+    }
+
+
+    // Determine neighbouring global cell Cells
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    labelListList neiGlobalCellCells(nBnd);
+    for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
+    {
+        neiGlobalCellCells[faceI-mesh.nInternalFaces()] =
+            globalCellCells[own[faceI]];
+    }
+    syncTools::swapBoundaryFaceList(mesh, neiGlobalCellCells, false);
+
+
+
+    // Construct stencil in global numbering
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    stencil_.setSize(mesh.nFaces());
+
+    labelHashSet faceStencil;
+
+    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+    {
+        faceStencil.clear();
+        label globalOwn = globalNumbering.toGlobal(own[faceI]);
+        faceStencil.insert(globalOwn);
+        const labelList& ownCCells = globalCellCells[own[faceI]];
+        forAll(ownCCells, i)
+        {
+            faceStencil.insert(ownCCells[i]);
+        }
+
+        label globalNei = globalNumbering.toGlobal(nei[faceI]);
+        faceStencil.insert(globalNei);
+        const labelList& neiCCells = globalCellCells[nei[faceI]];
+        forAll(neiCCells, i)
+        {
+            faceStencil.insert(neiCCells[i]);
+        }
+
+        // Guarantee owner first, neighbour second.
+        stencil_[faceI].setSize(faceStencil.size());
+        label n = 0;
+        stencil_[faceI][n++] = globalOwn;
+        stencil_[faceI][n++] = globalNei;
+        forAllConstIter(labelHashSet, faceStencil, iter)
+        {
+            if (iter.key() != globalOwn && iter.key() != globalNei)
+            {
+                stencil_[faceI][n++] = iter.key();
+            }
+        }
+        //Pout<< "internalface:" << faceI << " toc:" << faceStencil.toc()
+        //    << " stencil:" << stencil_[faceI] << endl;
+    }
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+        label faceI = pp.start();
+
+        if (pp.coupled())
+        {
+            forAll(pp, i)
+            {
+                faceStencil.clear();
+                label globalOwn = globalNumbering.toGlobal(own[faceI]);
+                faceStencil.insert(globalOwn);
+                const labelList& ownCCells = globalCellCells[own[faceI]];
+                forAll(ownCCells, i)
+                {
+                    faceStencil.insert(ownCCells[i]);
+                }
+                // Get the coupled cell
+                label globalNei = neiGlobal[faceI-mesh.nInternalFaces()];
+                faceStencil.insert(globalNei);
+                // And the neighbours of the coupled cell
+                const labelList& neiCCells =
+                    neiGlobalCellCells[faceI-mesh.nInternalFaces()];
+                forAll(neiCCells, i)
+                {
+                    faceStencil.insert(neiCCells[i]);
+                }
+
+                // Guarantee owner first, neighbour second.
+                stencil_[faceI].setSize(faceStencil.size());
+                label n = 0;
+                stencil_[faceI][n++] = globalOwn;
+                stencil_[faceI][n++] = globalNei;
+                forAllConstIter(labelHashSet, faceStencil, iter)
+                {
+                    if (iter.key() != globalOwn && iter.key() != globalNei)
+                    {
+                        stencil_[faceI][n++] = iter.key();
+                    }
+                }
+
+                //Pout<< "coupledface:" << faceI
+                //    << " toc:" << faceStencil.toc()
+                //    << " stencil:" << stencil_[faceI] << endl;
+
+                faceI++;
+            }
+        }
+        else if (!isA<emptyPolyPatch>(pp))
+        {
+            forAll(pp, i)
+            {
+                faceStencil.clear();
+                label globalOwn = globalNumbering.toGlobal(own[faceI]);
+                faceStencil.insert(globalOwn);
+                const labelList& ownCCells = globalCellCells[own[faceI]];
+                forAll(ownCCells, i)
+                {
+                    faceStencil.insert(ownCCells[i]);
+                }
+
+
+                // Guarantee owner first, neighbour second.
+                stencil_[faceI].setSize(faceStencil.size());
+                label n = 0;
+                stencil_[faceI][n++] = globalOwn;
+                forAllConstIter(labelHashSet, faceStencil, iter)
+                {
+                    if (iter.key() != globalOwn)
+                    {
+                        stencil_[faceI][n++] = iter.key();
+                    }
+                }
+
+                //Pout<< "boundaryface:" << faceI
+                //    << " toc:" << faceStencil.toc()
+                //    << " stencil:" << stencil_[faceI] << endl;
+
+                faceI++;
+            }
+        }
+    }
+}
+
+
+// Calculates extended stencil. This is per face
+// - owner
+// - cellCells of owner
+// - neighbour
+// - cellCells of neighbour
+// It comes in two parts:
+// - a map which collects/distributes all necessary data in a compact array
+// - the stencil (a labelList per face) which is a set of indices into this
+//   compact array.
+// The compact array is laid out as follows:
+// - first data for current processor (Pstream::myProcNo())
+//      - all cells
+//      - all boundary faces
+// - then per processor
+//      - all used cells and boundary faces
+void Foam::extendedStencil::calcExtendedFaceStencil(const polyMesh& mesh)
+{
+    const label nBnd = mesh.nFaces()-mesh.nInternalFaces();
+
+    // Global numbering for cells and boundary faces
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    globalIndex globalNumbering(mesh.nCells()+nBnd);
+
+
+    // Calculate stencil in global cell indices
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    calcFaceStencils(mesh, globalNumbering);
+
+
+    // Convert stencil to schedule
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // We now know what information we need from other processors. This needs
+    // to be converted into what information I need to send as well
+    // (mapDistribute)
+
+
+    // 1. Construct per processor compact addressing of the global cells
+    //    needed. The ones from the local processor are not included since
+    //    these are always all needed.
+    List<Map<label> > globalToProc(Pstream::nProcs());
+    {
+        const labelList& procPatchMap = mesh.globalData().procPatchMap();
+        const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+        // Presize with (as estimate) size of patch to neighbour.
+        forAll(procPatchMap, procI)
+        {
+            if (procPatchMap[procI] != -1)
+            {
+                globalToProc[procI].resize
+                (
+                    patches[procPatchMap[procI]].size()
+                );
+            }
+        }
+
+        // Collect all (non-local) globalcells/faces needed.
+        forAll(stencil_, faceI)
+        {
+            const labelList& stencilCells = stencil_[faceI];
+
+            forAll(stencilCells, i)
+            {
+                label globalCellI = stencilCells[i];
+                label procI = globalNumbering.whichProcID(stencilCells[i]);
+
+                if (procI != Pstream::myProcNo())
+                {
+                    label nCompact = globalToProc[procI].size();
+                    globalToProc[procI].insert(globalCellI, nCompact);
+                }
+            }
+        }
+        // Sort global cells needed (not really necessary)
+        forAll(globalToProc, procI)
+        {
+            if (procI != Pstream::myProcNo())
+            {
+                Map<label>& globalMap = globalToProc[procI];
+
+                SortableList<label> sorted(globalMap.toc());
+
+                forAll(sorted, i)
+                {
+                    Map<label>::iterator iter = globalMap.find(sorted[i]);
+                    iter() = i;
+                }
+            }
+        }
+
+
+        // forAll(globalToProc, procI)
+        // {
+        //     Pout<< "From processor:" << procI << " want cells/faces:" << endl;
+        //     forAllConstIter(Map<label>, globalToProc[procI], iter)
+        //     {
+        //         Pout<< "    global:" << iter.key()
+        //             << " local:" << globalNumbering.toLocal(procI, iter.key())
+        //             << endl;
+        //     }
+        //     Pout<< endl;
+        // }
+    }
+
+
+    // 2. The overall compact addressing is
+    // - myProcNo first
+    // - all other processors consecutively
+
+    labelList compactStart(Pstream::nProcs());
+    compactStart[Pstream::myProcNo()] = 0;
+    label nCompact = mesh.nCells()+nBnd;
+    forAll(compactStart, procI)
+    {
+        if (procI != Pstream::myProcNo())
+        {
+            compactStart[procI] = nCompact;
+            nCompact += globalToProc[procI].size();
+
+            // Pout<< "Data wanted from " << procI << " starts at "
+            //     << compactStart[procI] << endl;
+        }
+    }
+    // Pout<< "Overall cells needed:" << nCompact << endl;
+
+
+    // 3. Find out what to receive/send in compact addressing.
+    labelListList recvCompact(Pstream::nProcs());
+    for (label procI = 0; procI < Pstream::nProcs(); procI++)
+    {
+        if (procI != Pstream::myProcNo())
+        {
+            labelList wantedGlobals(globalToProc[procI].size());
+            recvCompact[procI].setSize(globalToProc[procI].size());
+
+            label i = 0;
+            forAllConstIter(Map<label>, globalToProc[procI], iter)
+            {
+                wantedGlobals[i] = iter.key();
+                recvCompact[procI][i] = compactStart[procI]+iter();
+                i++;
+            }
+
+            // Pout<< "From proc:" << procI
+            //     << " I need (globalcells):" << wantedGlobals
+            //     << " which are my compact:" << recvCompact[procI]
+            //     << endl;
+
+            // Send the global cell numbers I need from procI
+            OPstream str(Pstream::blocking, procI);
+            str << wantedGlobals;
+        }
+        else
+        {
+            recvCompact[procI] =
+                compactStart[procI]
+              + identity(mesh.nCells()+nBnd);
+        }
+    }
+    labelListList sendCompact(Pstream::nProcs());
+    for (label procI = 0; procI < Pstream::nProcs(); procI++)
+    {
+        if (procI != Pstream::myProcNo())
+        {
+            // See what neighbour wants to receive (= what I need to send)
+
+            IPstream str(Pstream::blocking, procI);
+            labelList globalCells(str);
+
+            labelList& procCompact = sendCompact[procI];
+            procCompact.setSize(globalCells.size());
+
+            // Convert from globalCells (all on my processor!) into compact
+            // addressing
+            forAll(globalCells, i)
+            {
+                label cellI = globalNumbering.toLocal(globalCells[i]);
+                procCompact[i] = compactStart[Pstream::myProcNo()]+cellI;
+            }
+        }
+        else
+        {
+            sendCompact[procI] = recvCompact[procI];
+        }
+    }
+
+    // Convert stencil to compact numbering
+    forAll(stencil_, faceI)
+    {
+        labelList& stencilCells = stencil_[faceI];
+
+        forAll(stencilCells, i)
+        {
+            label globalCellI = stencilCells[i];
+            label procI = globalNumbering.whichProcID(globalCellI);
+            if (procI != Pstream::myProcNo())
+            {
+                label localCompact = globalToProc[procI][globalCellI];
+                stencilCells[i] = compactStart[procI]+localCompact;
+            }
+            else
+            {
+                label localCompact = globalNumbering.toLocal(globalCellI);
+                stencilCells[i] = compactStart[procI]+localCompact;
+            }
+
+        }
+    }
+    // Pout<< "***stencil_:" << stencil_ << endl;
+
+    // Constuct map for distribution of compact data.
+    mapPtr_.reset
+    (
+        new mapDistribute
+        (
+            nCompact,
+            sendCompact,
+            recvCompact,
+            true            // reuse send/recv maps.
+        )
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::extendedStencil::extendedStencil
+(
+    const mapDistribute& map,
+    const labelListList& stencil
+)
+:
+    mapPtr_
+    (
+        autoPtr<mapDistribute>
+        (
+            new mapDistribute
+            (
+                map.constructSize(),
+                map.subMap(),
+                map.constructMap()
+            )
+        )
+    ),
+    stencil_(stencil)
+{}
+
+
+Foam::extendedStencil::extendedStencil(const polyMesh& mesh)
+{
+    calcExtendedFaceStencil(mesh);
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/extendedStencil.H

@@ -0,0 +1,151 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::extendedStencil
+
+Description
+    Calculates/constains the extended face stencil.
+
+    The stencil is a list of indices into either cells or boundary faces
+    in a compact way. (element 0 is owner, 1 is neighbour). The index numbering
+    is
+    - cells first
+    - then all (non-empty patch) boundary faces
+
+    When used in evaluation is a two stage process:
+    - collect the data (cell data and non-empty boundaries) into a
+    single field
+    - (parallel) distribute the field
+    - sum the weights*field.
+
+SourceFiles
+    extendedStencil.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef extendedStencil_H
+#define extendedStencil_H
+
+#include "mapDistribute.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class globalIndex;
+
+/*---------------------------------------------------------------------------*\
+                           Class extendedStencil Declaration
+\*---------------------------------------------------------------------------*/
+
+class extendedStencil
+{
+    // Private data
+
+        //- Swap map for getting neigbouring data
+        autoPtr<mapDistribute> mapPtr_;
+
+        //- Per face the stencil.
+        labelListList stencil_;
+
+
+    // Private Member Functions
+
+        void calcFaceStencils(const polyMesh&, const globalIndex&);
+
+        //- Calculate the stencil (but not weights)
+        void calcExtendedFaceStencil(const polyMesh&);
+
+
+        //- Disallow default bitwise copy construct
+        extendedStencil(const extendedStencil&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const extendedStencil&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        extendedStencil(const mapDistribute& map, const labelListList&);
+
+        //- Construct from all cells and boundary faces
+        extendedStencil(const polyMesh&);
+
+
+
+    // Member Functions
+
+        //- Return reference to the parallel distribution map
+        const mapDistribute& map() const
+        {
+            return mapPtr_();
+        }
+
+        //- Return reference to the stencil
+        const labelListList& stencil() const
+        {
+            return stencil_;
+        }
+
+        //- Use map to get the data into stencil order
+        template<class T>
+        void collectData
+        (
+            const GeometricField<T, fvPatchField, volMesh>& fld,
+            List<List<T> >& stencilFld
+        ) const;
+
+        //- Given weights interpolate vol field
+        template<class Type>
+        tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > interpolate
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& fld,
+            const List<List<scalar> >& stencilWeights
+        ) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "extendedStencilTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/extendedStencilTemplates.C

@@ -0,0 +1,129 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "extendedStencil.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::extendedStencil::collectData
+(
+    const GeometricField<Type, fvPatchField, volMesh>& fld,
+    List<List<Type> >& stencilFld
+) const
+{
+    // 1. Construct cell data in compact addressing
+    List<Type> compactFld(map().constructSize(), pTraits<Type>::zero);
+
+    // Insert my internal values
+    forAll(fld, cellI)
+    {
+        compactFld[cellI] = fld[cellI];
+    }
+    // Insert my boundary values
+    label nCompact = fld.size();
+    forAll(fld.boundaryField(), patchI)
+    {
+        const fvPatchField<Type>& pfld = fld.boundaryField()[patchI];
+
+        forAll(pfld, i)
+        {
+            compactFld[nCompact++] = pfld[i];
+        }
+    }
+
+    // Do all swapping
+    map().distribute(compactFld);
+
+    // 2. Pull to stencil
+    stencilFld.setSize(stencil_.size());
+
+    forAll(stencil_, faceI)
+    {
+        const labelList& compactCells = stencil_[faceI];
+
+        stencilFld[faceI].setSize(compactCells.size());
+
+        forAll(compactCells, i)
+        {
+            stencilFld[faceI][i] = compactFld[compactCells[i]];
+        }
+    }
+}
+
+
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> >
+Foam::extendedStencil::interpolate
+(
+    const GeometricField<Type, fvPatchField, volMesh>& fld,
+    const List<List<scalar> >& stencilWeights
+) const
+{
+    const fvMesh& mesh = fld.mesh();
+
+    // Collect internal and boundary values
+    List<List<Type> > stencilFld;
+    collectData(fld, stencilFld);
+
+    tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tsfCorr
+    (
+        new GeometricField<Type, fvsPatchField, surfaceMesh>
+        (
+            IOobject
+            (
+                fld.name(),
+                mesh.time().timeName(),
+                mesh
+            ),
+            mesh,
+            dimensioned<Type>
+            (
+                fld.name(),
+                fld.dimensions(),
+                pTraits<Type>::zero
+            )
+        )
+    );
+    GeometricField<Type, fvsPatchField, surfaceMesh>& sf = tsfCorr();
+
+    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+    {
+        const List<Type>& stField = stencilFld[faceI];
+        const List<scalar>& stWeight = stencilWeights[faceI];
+
+        forAll(stField, i)
+        {
+            sf[faceI] += stField[i]*stWeight[i];
+        }
+    }
+    // And what for boundaries?
+
+    return tsfCorr;
+}
+
+
+// ************************************************************************* //

src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.C

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

src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.H

@@ -0,0 +1,138 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    quadraticFit
+
+Description
+    Quadratic fit interpolation scheme which applies an explicit correction to
+    linear.
+
+SourceFiles
+    quadraticFit.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef quadraticFit_H
+#define quadraticFit_H
+
+#include "linear.H"
+#include "quadraticFitData.H"
+#include "extendedStencil.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class quadraticFit Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class quadraticFit
+:
+    public linear<Type>
+{
+    // Private Data
+        const scalar centralWeight_;
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        quadraticFit(const quadraticFit&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const quadraticFit&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("quadraticFit");
+
+
+    // Constructors
+
+        //- Construct from mesh and Istream
+        quadraticFit(const fvMesh& mesh, Istream& is)
+        :
+            linear<Type>(mesh),
+            centralWeight_(readScalar(is))
+        {}
+
+
+        //- Construct from mesh, faceFlux and Istream
+        quadraticFit
+        (
+            const fvMesh& mesh,
+            const surfaceScalarField& faceFlux,
+            Istream& is
+        )
+        :
+            linear<Type>(mesh),
+            centralWeight_(readScalar(is))
+        {}
+
+
+    // Member Functions
+
+        //- Return true if this scheme uses an explicit correction
+        virtual bool corrected() const
+        {
+            return true;
+        }
+
+        //- Return the explicit correction to the face-interpolate
+        virtual tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
+        correction
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& vf
+        ) const
+        {
+            const fvMesh& mesh = this->mesh();
+
+            const quadraticFitData& cfd = quadraticFitData::New
+            (
+                mesh,
+                centralWeight_
+            );
+
+            const extendedStencil& stencil = cfd.stencil();
+            const List<scalarList>& f = cfd.fit();
+
+            return stencil.interpolate(vf, f);
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C

@@ -0,0 +1,310 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "quadraticFitData.H"
+#include "surfaceFields.H"
+#include "volFields.H"
+#include "SVD.H"
+#include "syncTools.H"
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(quadraticFitData, 0);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::quadraticFitData::quadraticFitData
+(
+    const fvMesh& mesh,
+    const scalar cWeight
+)
+:
+    MeshObject<fvMesh, quadraticFitData>(mesh),
+    centralWeight_(cWeight),
+#   ifdef SPHERICAL_GEOMETRY
+    dim_(2),
+#   else
+    dim_(mesh.nGeometricD()),
+#   endif
+    minSize_
+    (
+        dim_ == 1 ? 3 :
+        dim_ == 2 ? 6 :
+        dim_ == 3 ? 9 : 0
+    ),
+    stencil_(mesh),
+    fit_(mesh.nInternalFaces())
+{
+    if (debug)
+    {
+        Info << "Contructing quadraticFitData" << endl;
+    }
+
+    // check input
+    if (centralWeight_ < 1 - SMALL)
+    {
+        FatalErrorIn("quadraticFitData::quadraticFitData")
+            << "centralWeight requested = " << centralWeight_
+            << " should not be less than one"
+            << exit(FatalError);
+    }
+
+    if (minSize_ == 0)
+    {
+        FatalErrorIn("quadraticFitSnGradData")
+            << " dimension must be 1,2 or 3, not" << dim_ << exit(FatalError);
+    }
+
+    // store the polynomial size for each cell to write out
+    surfaceScalarField interpPolySize
+    (
+        IOobject
+        (
+            "quadraticFitInterpPolySize",
+            "constant",
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("quadraticFitInterpPolySize", dimless, scalar(0))
+    );
+
+    // Get the cell/face centres in stencil order.
+    // Centred face stencils no good for triangles of tets. Need bigger stencils
+    List<List<point> > stencilPoints(stencil_.stencil().size());
+    stencil_.collectData
+    (
+        mesh.C(),
+        stencilPoints
+    );
+
+    // find the fit coefficients for every face in the mesh
+
+    for(label faci = 0; faci < mesh.nInternalFaces(); faci++)
+    {
+        interpPolySize[faci] = calcFit(stencilPoints[faci], faci);
+    }
+    interpPolySize.write();
+    if (debug)
+    {
+        Info<< "quadraticFitData::quadraticFitData() :"
+            << "Finished constructing polynomialFit data"
+            << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::quadraticFitData::findFaceDirs
+(
+    vector& idir,        // value changed in return
+    vector& jdir,        // value changed in return
+    vector& kdir,        // value changed in return
+    const fvMesh& mesh,
+    const label faci
+)
+{
+    idir = mesh.Sf()[faci];
+    idir /= mag(idir);
+
+#   ifndef SPHERICAL_GEOMETRY
+    if (mesh.nGeometricD() <= 2) // find the normal direcion
+    {
+        if (mesh.directions()[0] == -1)
+        {
+            kdir = vector(1, 0, 0);
+        }
+        else if (mesh.directions()[1] == -1)
+        {
+            kdir = vector(0, 1, 0);
+        }
+        else
+        {
+            kdir = vector(0, 0, 1);
+        }
+    }
+    else // 3D so find a direction in the place of the face
+    {
+        const face& f = mesh.faces()[faci];
+        kdir = mesh.points()[f[0]] - mesh.points()[f[1]];
+    }
+#   else
+    // Spherical geometry so kdir is the radial direction
+    kdir = mesh.Cf()[faci];
+#   endif
+
+    if (mesh.nGeometricD() == 3)
+    {
+        // Remove the idir component from kdir and normalise
+        kdir -= (idir & kdir)*idir;
+
+        scalar magk = mag(kdir);
+
+        if (magk < SMALL)
+        {
+            FatalErrorIn("findFaceDirs") << " calculated kdir = zero"
+                << exit(FatalError);
+        }
+        else
+        {
+            kdir /= magk;
+        }
+    }
+
+    jdir = kdir ^ idir;
+}
+
+
+Foam::label Foam::quadraticFitData::calcFit
+(
+    const List<point>& C,
+    const label faci
+)
+{
+    vector idir(1,0,0);
+    vector jdir(0,1,0);
+    vector kdir(0,0,1);
+    findFaceDirs(idir, jdir, kdir, mesh(), faci);
+
+    scalarList wts(C.size(), scalar(1));
+    wts[0] = centralWeight_;
+    wts[1] = centralWeight_;
+
+    point p0 = mesh().faceCentres()[faci];
+    scalar scale = 0;
+
+    // calculate the matrix of the polynomial components
+    scalarRectangularMatrix B(C.size(), minSize_, scalar(0));
+
+    for(label ip = 0; ip < C.size(); ip++)
+    {
+        const point& p = C[ip];
+
+        scalar px = (p - p0)&idir;
+        scalar py = (p - p0)&jdir;
+#       ifndef SPHERICAL_GEOMETRY
+        scalar pz = (p - p0)&kdir;
+#       else
+        scalar pz = mag(p) - mag(p0);
+#       endif
+
+        if (ip == 0)
+        {
+            scale = max(max(mag(px), mag(py)), mag(pz));
+        }
+
+        px /= scale;
+        py /= scale;
+        pz /= scale;
+
+        label is = 0;
+        B[ip][is++] = wts[ip];
+
+        B[ip][is++] = wts[ip]*px;
+        B[ip][is++] = wts[ip]*sqr(px);
+
+        if (dim_ >= 2)
+        {
+            B[ip][is++] = wts[ip]*py;
+            B[ip][is++] = wts[ip]*px*py;
+            B[ip][is++] = wts[ip]*sqr(py);
+        }
+        if (dim_ == 3)
+        {
+            B[ip][is++] = wts[ip]*pz;
+            B[ip][is++] = wts[ip]*px*pz;
+            //B[ip][is++] = wts[ip]*py*pz;
+            B[ip][is++] = wts[ip]*sqr(pz);
+        }
+    }
+
+    // Set the fit
+    label stencilSize = C.size();
+    fit_[faci].setSize(stencilSize);
+    scalarList singVals(minSize_);
+    label nSVDzeros = 0;
+
+    bool goodFit = false;
+    for(scalar tol = SMALL; tol < 0.1 && !goodFit; tol *= 10)
+    {
+        SVD svd(B, tol);
+
+        scalar fit0 = wts[0]*svd.VSinvUt()[0][0];
+        scalar fit1 = wts[1]*svd.VSinvUt()[0][1];
+
+        goodFit = sign(fit0) == sign(fit1);
+
+        if (goodFit)
+        {
+            fit_[faci][0] = fit0;
+            fit_[faci][1] = fit1;
+            for(label i = 2; i < stencilSize; i++)
+            {
+                fit_[faci][i] = wts[i]*svd.VSinvUt()[0][i];
+            }
+
+            singVals = svd.S();
+            nSVDzeros = svd.nZeros();
+        }
+    }
+    if (!goodFit)
+    {
+        FatalErrorIn
+        (
+            "quadraticFitData::calcFit(const pointField&, const label)"
+            ) << "For face " << faci << endl
+            << "Fit not good even once tol >= 0.1"
+            << exit(FatalError);
+    }
+
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>& w =
+        mesh().surfaceInterpolation::weights();
+
+    // remove the uncorrected linear coefficients
+
+    fit_[faci][0] -= w[faci];
+    fit_[faci][1] -= 1 - w[faci];
+
+    return minSize_ - nSVDzeros;
+}
+
+
+bool Foam::quadraticFitData::movePoints()
+{
+    notImplemented("quadraticFitData::movePoints()");
+
+    return true;
+}
+
+
+// ************************************************************************* //

src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.H

@@ -0,0 +1,140 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    quadraticFitData
+
+Description
+    Data for the quadratic fit correction interpolation scheme
+
+SourceFiles
+    quadraticFitData.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef quadraticFitData_H
+#define quadraticFitData_H
+
+#include "MeshObject.H"
+#include "fvMesh.H"
+#include "extendedStencil.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class globalIndex;
+
+/*---------------------------------------------------------------------------*\
+                    Class quadraticFitData Declaration
+\*---------------------------------------------------------------------------*/
+
+class quadraticFitData
+:
+    public MeshObject<fvMesh, quadraticFitData>
+{
+    // Private data
+
+        //- weights for central stencil
+        const scalar centralWeight_;
+
+        //- dimensionality of the geometry
+        const label dim_;
+
+        //- minimum stencil size
+        const label minSize_;
+
+        //- Extended stencil addressing
+        extendedStencil stencil_;
+
+        //- For each cell in the mesh store the values which multiply the
+        //  values of the stencil to obtain the gradient for each direction
+        List<scalarList> fit_;
+
+
+    // Private member functions
+
+        //- Find the normal direction and i, j and k directions for face faci
+        static void findFaceDirs
+        (
+            vector& idir,        // value changed in return
+            vector& jdir,        // value changed in return
+            vector& kdir,        // value changed in return
+            const fvMesh& mesh,
+            const label faci
+        );
+
+        label calcFit(const List<point>&, const label faci);
+
+
+public:
+
+    TypeName("quadraticFitData");
+
+
+    // Constructors
+
+        explicit quadraticFitData
+        (
+            const fvMesh& mesh,
+            scalar cWeightDim
+        );
+
+
+    // Destructor
+
+        virtual ~quadraticFitData()
+        {}
+
+
+    // Member functions
+
+
+        //- Return reference to the stencil
+        const extendedStencil& stencil() const
+        {
+            return stencil_;
+        }
+
+        //- Return reference to fit coefficients
+        const List<scalarList>& fit() const
+        {
+            return fit_;
+        }
+
+        //- Delete the data when the mesh moves not implemented
+        virtual bool movePoints();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.C

@@ -114,7 +114,7 @@ Foam::scalarField Foam::chemistryModel::omega
     forAll(reactions_, i)
     {
         const reaction& R = reactions_[i];
-        
+
         scalar omegai = omega
         (
             R, c, T, p, pf, cf, lRef, pr, cr, rRef
@@ -164,13 +164,13 @@ Foam::scalar Foam::chemistryModel::omega
 
     pf = 1.0;
     pr = 1.0;
-    
+
     label Nl = R.lhs().size();
     label Nr = R.rhs().size();
-    
+
     label slRef = 0;
     lRef = R.lhs()[slRef].index;
-        
+
     pf = kf;
     for(label s=1; s<Nl; s++)
     {
@@ -212,7 +212,7 @@ Foam::scalar Foam::chemistryModel::omega
 
     label srRef = 0;
     rRef = R.rhs()[srRef].index;
-    
+
     // find the matrix element and element position for the rhs
     pr = kr;
     for(label s=1; s<Nr; s++)
@@ -250,7 +250,7 @@ Foam::scalar Foam::chemistryModel::omega
         {
             pr *= pow(cr, exp-1.0);
         }
-        
+
     }
 
     return pf*cf - pr*cr;
@@ -313,12 +313,12 @@ void Foam::chemistryModel::jacobian
     const scalar t,
     const scalarField& c,
     scalarField& dcdt,
-    Matrix<scalar>& dfdc
+    scalarSquareMatrix& dfdc
 ) const
 {
     scalar T = c[Ns_];
     scalar p = c[Ns_ + 1];
-    
+
     scalarField c2(Ns(), 0.0);
     for(label i=0; i<Ns(); i++)
     {
@@ -470,23 +470,23 @@ Foam::tmp<Foam::volScalarField> Foam::chemistryModel::tc() const
             scalar pi = thermo_.p()[celli];
             scalarField c(Ns_);
             scalar cSum = 0.0;
-            
+
             for(label i=0; i<Ns_; i++)
             {
                 scalar Yi = Y_[i][celli];
                 c[i] = rhoi*Yi/specieThermo_[i].W();
                 cSum += c[i];
             }
-            
+
             forAll(reactions_, i)
             {
                 const reaction& R = reactions_[i];
-                
+
                 omega
                 (
                     R, c, Ti, pi, pf, cf, lRef, pr, cr, rRef
                 );
-                
+
                 forAll(R.rhs(), s)
                 {
                     scalar sr = R.rhs()[s].stoichCoeff;
@@ -544,22 +544,22 @@ void Foam::chemistryModel::calculate()
             {
                 RR_[i][celli] = 0.0;
             }
-            
+
             scalar rhoi = rho_[celli];
             scalar Ti = thermo_.T()[celli];
             scalar pi = thermo_.p()[celli];
-            
+
             scalarField c(Ns_);
             scalarField dcdt(nEqns(), 0.0);
-            
+
             for(label i=0; i<Ns_; i++)
             {
                 scalar Yi = Y_[i][celli];
                 c[i] = rhoi*Yi/specieThermo_[i].W();
             }
-            
+
             dcdt = omega(c, Ti, pi);
-            
+
             for(label i=0; i<Ns_; i++)
             {
                 RR_[i][celli] = dcdt[i]*specieThermo_[i].W();
@@ -624,7 +624,7 @@ Foam::scalar Foam::chemistryModel::solve(const scalar t0, const scalar deltaT)
             for(label i=0; i<Ns_; i++)
             {
                 mixture += (c[i]/cTot)*specieThermo_[i];
-            }        
+            }
             Ti = mixture.TH(hi, Ti);
 
             timeLeft -= dt;
@@ -639,7 +639,7 @@ Foam::scalar Foam::chemistryModel::solve(const scalar t0, const scalar deltaT)
         for(label i=0; i<Ns_; i++)
         {
             WTot += c[i]*specieThermo_[i].W();
-        }        
+        }
         WTot /= cTot;
 
         for(label i=0; i<Ns_; i++)

src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.H

@@ -39,7 +39,6 @@ SourceFiles
 
 #include "hCombustionThermo.H"
 #include "reactingMixture.H"
-#include "Matrix.H"
 #include "ODE.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -251,7 +250,7 @@ public:
             const scalar t,
             const scalarField& c,
             scalarField& dcdt,
-            Matrix<scalar>& dfdc
+            scalarSquareMatrix& dfdc
         ) const;
 
         //- Calculates the reaction rates

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

@@ -204,11 +204,20 @@ LRR::LRR
             IOobject::AUTO_WRITE
         ),
         autoCreateMut("mut", mesh_)
+    ),
+    alphat_
+    (
+        IOobject
+        (
+            "alphat",
+            runTime_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        autoCreateAlphat("alphat", mesh_)
     )
 {
-    mut_ = Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
-    mut_.correctBoundaryConditions();
-
     if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
     {
         FatalErrorIn
@@ -221,6 +230,12 @@ LRR::LRR
             << exit(FatalError);
     }
 
+    mut_ == Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    mut_.correctBoundaryConditions();
+
+    alphat_ == mut_/Prt_;
+    alphat_.correctBoundaryConditions();
+
     printCoeffs();
 }
 
@@ -310,8 +325,13 @@ void LRR::correct()
     if (!turbulence_)
     {
         // Re-calculate viscosity
-        mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
+        mut_ == rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
         mut_.correctBoundaryConditions();
+
+        // Re-calculate thermal diffusivity
+        alphat_ = mut_/Prt_;
+        alphat_.correctBoundaryConditions();
+
         return;
     }
 
@@ -399,9 +419,13 @@ void LRR::correct()
 
 
     // Re-calculate viscosity
-    mut_ = rho_*Cmu_*sqr(k_)/epsilon_;
+    mut_ == rho_*Cmu_*sqr(k_)/epsilon_;
     mut_.correctBoundaryConditions();
 
+    // Re-calculate thermal diffusivity
+    alphat_ = mut_/Prt_;
+    alphat_.correctBoundaryConditions();
+
 
     // Correct wall shear stresses
 

src/turbulenceModels/RAS/compressible/LRR/LRR.H

@@ -97,6 +97,7 @@ class LRR
         volScalarField k_;
         volScalarField epsilon_;
         volScalarField mut_;
+        volScalarField alphat_;
 
 
 public:
@@ -152,7 +153,7 @@ public:
         {
             return tmp<volScalarField>
             (
-                new volScalarField("alphaEff", alphah_*mut_ + alpha())
+                new volScalarField("alphaEff", alphah_*alphat_ + alpha())
             );
         }
 

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

@@ -226,11 +226,20 @@ LaunderGibsonRSTM::LaunderGibsonRSTM
             IOobject::AUTO_WRITE
         ),
         autoCreateMut("mut", mesh_)
+    ),
+    alphat_
+    (
+        IOobject
+        (
+            "alphat",
+            runTime_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        autoCreateAlphat("alphat", mesh_)
     )
 {
-    mut_ = Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
-    mut_.correctBoundaryConditions();
-
     if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
     {
         FatalErrorIn
@@ -243,6 +252,12 @@ LaunderGibsonRSTM::LaunderGibsonRSTM
             << exit(FatalError);
     }
 
+    mut_ == Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    mut_.correctBoundaryConditions();
+
+    alphat_ == mut_/Prt_;
+    alphat_.correctBoundaryConditions();
+
     printCoeffs();
 }
 
@@ -335,8 +350,13 @@ void LaunderGibsonRSTM::correct()
     if (!turbulence_)
     {
         // Re-calculate viscosity
-        mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
+        mut_ == rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
         mut_.correctBoundaryConditions();
+
+        // Re-calculate thermal diffusivity
+        alphat_ = mut_/Prt_;
+        alphat_.correctBoundaryConditions();
+
         return;
     }
 
@@ -438,9 +458,12 @@ void LaunderGibsonRSTM::correct()
 
 
     // Re-calculate turbulent viscosity
-    mut_ = Cmu_*rho_*sqr(k_)/epsilon_;
+    mut_ == Cmu_*rho_*sqr(k_)/epsilon_;
     mut_.correctBoundaryConditions();
 
+    // Re-calculate thermal diffusivity
+    alphat_ = mut_/Prt_;
+    alphat_.correctBoundaryConditions();
 
     // Correct wall shear stresses
 

src/turbulenceModels/RAS/compressible/LaunderGibsonRSTM/LaunderGibsonRSTM.H

@@ -104,6 +104,7 @@ class LaunderGibsonRSTM
         volScalarField k_;
         volScalarField epsilon_;
         volScalarField mut_;
+        volScalarField alphat_;
 
 
 public:
@@ -161,7 +162,7 @@ public:
         {
             return tmp<volScalarField>
             (
-                new volScalarField("alphaEff", alphah_*mut_ + alpha())
+                new volScalarField("alphaEff", alphah_*alphat_ + alpha())
             );
         }
 

src/turbulenceModels/RAS/compressible/Make/files

@@ -14,6 +14,9 @@ kOmegaSST/kOmegaSST.C
 /* Wall functions */
 wallFunctions = derivedFvPatchFields/wallFunctions
 
+alphatWallFunctions = $(wallFunctions)/alphatWallFunctions
+$(alphatWallFunctions)/alphatWallFunction/alphatWallFunctionFvPatchScalarField.C
+
 mutWallFunctions = $(wallFunctions)/mutWallFunctions
 $(mutWallFunctions)/mutWallFunction/mutWallFunctionFvPatchScalarField.C
 $(mutWallFunctions)/mutRoughWallFunction/mutRoughWallFunctionFvPatchScalarField.C

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

@@ -47,7 +47,8 @@ void RASModel::printCoeffs()
 {
     if (printCoeffs_)
     {
-        Info<< type() << "Coeffs" << coeffDict_ << endl;
+        Info<< type() << "Coeffs" << coeffDict_ << nl
+            << "wallFunctionCoeffs" << wallFunctionDict_ << endl;
     }
 }
 
@@ -115,6 +116,15 @@ RASModel::RASModel
             0.09
         )
     ),
+    Prt_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Prt",
+            wallFunctionDict_,
+            0.85
+        )
+    ),
 
     yPlusLam_(yPlusLam(kappa_.value(), E_.value())),
 
@@ -148,11 +158,9 @@ tmp<scalarField> RASModel::yPlus(const label patchNo) const
     tmp<scalarField> tYp(new scalarField(curPatch.size()));
     scalarField& Yp = tYp();
 
-    if (typeid(curPatch) == typeid(wallFvPatch))
+    if (isType<wallFvPatch>(curPatch))
     {
-        scalar Cmu(readScalar(coeffDict_.lookup("Cmu")));
-
-        Yp = pow(Cmu, 0.25)
+        Yp = pow(Cmu_.value(), 0.25)
             *y_[patchNo]
             *sqrt(k()().boundaryField()[patchNo].patchInternalField())
            /(
@@ -165,8 +173,8 @@ tmp<scalarField> RASModel::yPlus(const label patchNo) const
         WarningIn
         (
             "tmp<scalarField> RASModel::yPlus(const label patchNo) const"
-        )   << "Patch " << patchNo << " is not a wall.  Returning blank field"
-            << endl;
+        )   << "Patch " << patchNo << " is not a wall. Returning null field"
+            << nl << endl;
 
         Yp.setSize(0);
     }
@@ -191,8 +199,11 @@ bool RASModel::read()
         lookup("turbulence") >> turbulence_;
         coeffDict_ = subDict(type() + "Coeffs");
 
-        kappa_.readIfPresent(subDict("wallFunctionCoeffs"));
-        E_.readIfPresent(subDict("wallFunctionCoeffs"));
+        wallFunctionDict_ = subDict("wallFunctionCoeffs");
+        kappa_.readIfPresent(wallFunctionDict_);
+        E_.readIfPresent(wallFunctionDict_);
+        Cmu_.readIfPresent(wallFunctionDict_);
+        Prt_.readIfPresent(wallFunctionDict_);
 
         yPlusLam_ = yPlusLam(kappa_.value(), E_.value());
 

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

@@ -95,6 +95,7 @@ protected:
         dimensionedScalar kappa_;
         dimensionedScalar E_;
         dimensionedScalar Cmu_;
+        dimensionedScalar Prt_;
 
         scalar yPlusLam_;
 
@@ -244,6 +245,12 @@ public:
                 return Cmu_;
             }
 
+            //- Return turbulent Prandtl number for use in wall-functions
+            dimensionedScalar Prt() const
+            {
+                return Prt_;
+            }
+
             //- Return the near wall distances
             const nearWallDist& y() const
             {

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

@@ -174,11 +174,26 @@ RNGkEpsilon::RNGkEpsilon
             IOobject::AUTO_WRITE
         ),
         autoCreateMut("mut", mesh_)
+    ),
+    alphat_
+    (
+        IOobject
+        (
+            "alphat",
+            runTime_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        autoCreateAlphat("alphat", mesh_)
     )
 {
-    mut_ = Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    mut_ == Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
     mut_.correctBoundaryConditions();
 
+    alphat_ == mut_/Prt_;
+    alphat_.correctBoundaryConditions();
+
     printCoeffs();
 }
 
@@ -263,8 +278,13 @@ void RNGkEpsilon::correct()
     if (!turbulence_)
     {
         // Re-calculate viscosity
-        mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
+        mut_ == rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
         mut_.correctBoundaryConditions();
+
+        // Re-calculate thermal diffusivity
+        alphat_ = mut_/Prt_;
+        alphat_.correctBoundaryConditions();
+
         return;
     }
 
@@ -330,8 +350,12 @@ void RNGkEpsilon::correct()
 
 
     // Re-calculate viscosity
-    mut_ = rho_*Cmu_*sqr(k_)/epsilon_;
+    mut_ == rho_*Cmu_*sqr(k_)/epsilon_;
     mut_.correctBoundaryConditions();
+
+    // Re-calculate thermal diffusivity
+    alphat_ = mut_/Prt_;
+    alphat_.correctBoundaryConditions();
 }
 
 

src/turbulenceModels/RAS/compressible/RNGkEpsilon/RNGkEpsilon.H

@@ -87,6 +87,7 @@ class RNGkEpsilon
         volScalarField k_;
         volScalarField epsilon_;
         volScalarField mut_;
+        volScalarField alphat_;
 
 
 public:
@@ -142,7 +143,7 @@ public:
         {
             return tmp<volScalarField>
             (
-                new volScalarField("alphaEff", alphah_*mut_ + alpha())
+                new volScalarField("alphaEff", alphah_*alphat_ + alpha())
             );
         }
 

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

@@ -27,6 +27,7 @@ License
 #include "backwardsCompatibilityWallFunctions.H"
 
 #include "calculatedFvPatchField.H"
+#include "alphatWallFunctionFvPatchScalarField.H"
 #include "mutWallFunctionFvPatchScalarField.H"
 #include "epsilonWallFunctionFvPatchScalarField.H"
 #include "kQRWallFunctionFvPatchField.H"
@@ -41,6 +42,76 @@ namespace compressible
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+tmp<volScalarField> autoCreateAlphat
+(
+    const word& fieldName,
+    const fvMesh& mesh
+)
+{
+    IOobject alphatHeader
+    (
+        fieldName,
+        mesh.time().timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE,
+        false
+    );
+
+    if (alphatHeader.headerOk())
+    {
+        return tmp<volScalarField>(new volScalarField(alphatHeader, mesh));
+    }
+    else
+    {
+        Info<< "--> Upgrading " << fieldName << " to employ run-time "
+            << "selectable wall functions" << endl;
+
+        const fvBoundaryMesh& bm = mesh.boundary();
+
+        wordList alphatBoundaryTypes(bm.size());
+
+        forAll(bm, patchI)
+        {
+            if (isType<wallFvPatch>(bm[patchI]))
+            {
+                alphatBoundaryTypes[patchI] =
+                    RASModels::alphatWallFunctionFvPatchScalarField::typeName;
+            }
+            else
+            {
+                alphatBoundaryTypes[patchI] =
+                    calculatedFvPatchField<scalar>::typeName;
+            }
+        }
+
+        tmp<volScalarField> alphat
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    fieldName,
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                mesh,
+                dimensionedScalar("zero", dimDensity*dimArea/dimTime, 0.0),
+                alphatBoundaryTypes
+            )
+        );
+
+        Info<< "    Writing updated " << fieldName << endl;
+        alphat().write();
+
+        return alphat;
+    }
+}
+
+
 tmp<volScalarField> autoCreateMut
 (
     const word& fieldName,

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

@@ -53,6 +53,13 @@ namespace compressible
         const fvMesh& mesh
     );
 
+    //- alphat
+    tmp<volScalarField> autoCreateAlphat
+    (
+        const word& fieldName,
+        const fvMesh& mesh
+    );
+
     //- epsilon
     tmp<volScalarField> autoCreateEpsilon
     (

src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatWallFunction/alphatWallFunctionFvPatchScalarField.C

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

src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatWallFunction/alphatWallFunctionFvPatchScalarField.H

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

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

@@ -155,11 +155,26 @@ kEpsilon::kEpsilon
             IOobject::AUTO_WRITE
         ),
         autoCreateMut("mut", mesh_)
+    ),
+    alphat_
+    (
+        IOobject
+        (
+            "alphat",
+            runTime_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        autoCreateAlphat("alphat", mesh_)
     )
 {
-    mut_ = Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    mut_ == Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
     mut_.correctBoundaryConditions();
 
+    alphat_ == mut_/Prt_;
+    alphat_.correctBoundaryConditions();
+
     printCoeffs();
 }
 
@@ -243,8 +258,13 @@ void kEpsilon::correct()
     if (!turbulence_)
     {
         // Re-calculate viscosity
-        mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
+        mut_ == rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
         mut_.correctBoundaryConditions();
+
+        // Re-calculate thermal diffusivity
+        alphat_ = mut_/Prt_;
+        alphat_.correctBoundaryConditions();
+
         return;
     }
 
@@ -303,8 +323,12 @@ void kEpsilon::correct()
 
 
     // Re-calculate viscosity
-    mut_ = rho_*Cmu_*sqr(k_)/epsilon_;
+    mut_ == rho_*Cmu_*sqr(k_)/epsilon_;
     mut_.correctBoundaryConditions();
+
+    // Re-calculate thermal diffusivity
+    alphat_ = mut_/Prt_;
+    alphat_.correctBoundaryConditions();
 }
 
 

src/turbulenceModels/RAS/compressible/kEpsilon/kEpsilon.H

@@ -74,7 +74,6 @@ class kEpsilon
 
         // Model coefficients
 
-//            dimensionedScalar Cmu;
             dimensionedScalar Cmu_;
             dimensionedScalar C1_;
             dimensionedScalar C2_;
@@ -88,6 +87,7 @@ class kEpsilon
             volScalarField k_;
             volScalarField epsilon_;
             volScalarField mut_;
+            volScalarField alphat_;
 
 
 public:
@@ -144,7 +144,7 @@ public:
         {
             return tmp<volScalarField>
             (
-                new volScalarField("alphaEff", alphah_*mut_ + alpha())
+                new volScalarField("alphaEff", alphah_*alphat_ + alpha())
             );
         }
 

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

@@ -258,11 +258,26 @@ kOmegaSST::kOmegaSST
             IOobject::AUTO_WRITE
         ),
         autoCreateMut("mut", mesh_)
+    ),
+    alphat_
+    (
+        IOobject
+        (
+            "alphat",
+            runTime_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        autoCreateAlphat("alphat", mesh_)
     )
 {
-    mut_ = a1_*rho_*k_/max(a1_*omega_, F2()*sqrt(magSqr(symm(fvc::grad(U_)))));
+    mut_ == a1_*rho_*k_/max(a1_*omega_, F2()*sqrt(magSqr(symm(fvc::grad(U_)))));
     mut_.correctBoundaryConditions();
 
+    alphat_ == mut_/Prt_;
+    alphat_.correctBoundaryConditions();
+
     printCoeffs();
 }
 
@@ -351,11 +366,15 @@ void kOmegaSST::correct()
     if (!turbulence_)
     {
         // Re-calculate viscosity
-        mut_ =
+        mut_ ==
             a1_*rho_*k_
            /max(a1_*omega_, F2()*sqrt(magSqr(symm(fvc::grad(U_)))));
         mut_.correctBoundaryConditions();
 
+        // Re-calculate thermal diffusivity
+        alphat_ = mut_/Prt_;
+        alphat_.correctBoundaryConditions();
+
         return;
     }
 
@@ -430,8 +449,12 @@ void kOmegaSST::correct()
 
 
     // Re-calculate viscosity
-    mut_ = a1_*rho_*k_/max(a1_*omega_, F2()*sqrt(S2));
+    mut_ == a1_*rho_*k_/max(a1_*omega_, F2()*sqrt(S2));
     mut_.correctBoundaryConditions();
+
+    // Re-calculate thermal diffusivity
+    alphat_ = mut_/Prt_;
+    alphat_.correctBoundaryConditions();
 }
 
 

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

@@ -134,6 +134,7 @@ class kOmegaSST
         volScalarField k_;
         volScalarField omega_;
         volScalarField mut_;
+        volScalarField alphat_;
 
 
     // Private member functions
@@ -238,7 +239,7 @@ public:
         {
             return tmp<volScalarField>
             (
-                new volScalarField("alphaEff", alphah_*mut_ + alpha())
+                new volScalarField("alphaEff", alphah_*alphat_ + alpha())
             );
         }
 

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

@@ -187,14 +187,29 @@ realizableKE::realizableKE
             IOobject::AUTO_WRITE
         ),
         autoCreateMut("mut", mesh_)
+    ),
+    alphat_
+    (
+        IOobject
+        (
+            "alphat",
+            runTime_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        autoCreateAlphat("alphat", mesh_)
     )
 {
     bound(k_, k0_);
     bound(epsilon_, epsilon0_);
 
-    mut_ = rCmu(fvc::grad(U_))*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    mut_ == rCmu(fvc::grad(U_))*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
     mut_.correctBoundaryConditions();
 
+    alphat_ == mut_/Prt_;
+    alphat_.correctBoundaryConditions();
+
     printCoeffs();
 }
 
@@ -276,8 +291,13 @@ void realizableKE::correct()
     if (!turbulence_)
     {
         // Re-calculate viscosity
-        mut_ = rCmu(fvc::grad(U_))*rho_*sqr(k_)/epsilon_;
+        mut_ == rCmu(fvc::grad(U_))*rho_*sqr(k_)/epsilon_;
         mut_.correctBoundaryConditions();
+
+        // Re-calculate thermal diffusivity
+        alphat_ = mut_/Prt_;
+        alphat_.correctBoundaryConditions();
+
         return;
     }
 
@@ -342,8 +362,12 @@ void realizableKE::correct()
     bound(k_, k0_);
 
     // Re-calculate viscosity
-    mut_ = rCmu(gradU, S2, magS)*rho_*sqr(k_)/epsilon_;
+    mut_ == rCmu(gradU, S2, magS)*rho_*sqr(k_)/epsilon_;
     mut_.correctBoundaryConditions();
+
+    // Re-calculate thermal diffusivity
+    alphat_ = mut_/Prt_;
+    alphat_.correctBoundaryConditions();
 }
 
 

src/turbulenceModels/RAS/compressible/realizableKE/realizableKE.H

@@ -91,6 +91,7 @@ class realizableKE
         volScalarField k_;
         volScalarField epsilon_;
         volScalarField mut_;
+        volScalarField alphat_;
 
         tmp<volScalarField> rCmu
         (
@@ -157,7 +158,7 @@ public:
         {
             return tmp<volScalarField>
             (
-                new volScalarField("alphaEff", alphah_*mut_ + alpha())
+                new volScalarField("alphaEff", alphah_*alphat_ + alpha())
             );
         }
 

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

@@ -187,7 +187,7 @@ LRR::LRR
         autoCreateNut("nut", mesh_)
     )
 {
-    nut_ = Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    nut_ == Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
     nut_.correctBoundaryConditions();
 
     if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
@@ -381,7 +381,7 @@ void LRR::correct()
 
 
     // Re-calculate viscosity
-    nut_ = Cmu_*sqr(k_)/epsilon_;
+    nut_ == Cmu_*sqr(k_)/epsilon_;
     nut_.correctBoundaryConditions();
 
 

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

@@ -216,7 +216,7 @@ LaunderGibsonRSTM::LaunderGibsonRSTM
         autoCreateNut("nut", mesh_)
     )
 {
-    nut_ = Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    nut_ == Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
     nut_.correctBoundaryConditions();
 
     if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
@@ -422,7 +422,7 @@ void LaunderGibsonRSTM::correct()
 
 
     // Re-calculate turbulent viscosity
-    nut_ = Cmu_*sqr(k_)/epsilon_;
+    nut_ == Cmu_*sqr(k_)/epsilon_;
     nut_.correctBoundaryConditions();
 
 

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

@@ -228,7 +228,7 @@ LienCubicKE::LienCubicKE
         )
     )
 {
-    nut_ = Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_) + C5viscosity_;
+    nut_ == Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_) + C5viscosity_;
     nut_.correctBoundaryConditions();
 
     printCoeffs();
@@ -385,7 +385,7 @@ void LienCubicKE::correct()
         - 2.0*pow(Cmu_, 3.0)*pow(k_, 4.0)/pow(epsilon_, 3.0)
        *(magSqr(gradU_ + gradU_.T()) - magSqr(gradU_ - gradU_.T()));
 
-    nut_ = Cmu_*sqr(k_)/epsilon_ + C5viscosity_;
+    nut_ == Cmu_*sqr(k_)/epsilon_ + C5viscosity_;
     nut_.correctBoundaryConditions();
 
     nonlinearStress_ = symm

src/turbulenceModels/RAS/incompressible/RASModel/RASModel.C

@@ -45,7 +45,8 @@ void RASModel::printCoeffs()
 {
     if (printCoeffs_)
     {
-        Info<< type() << "Coeffs" << coeffDict_ << endl;;
+        Info<< type() << "Coeffs" << coeffDict_ << nl
+            << "wallFunctionCoeffs" << wallFunctionDict_ << endl;
     }
 }
 
@@ -148,9 +149,10 @@ tmp<scalarField> RASModel::yPlus(const label patchNo) const
     tmp<scalarField> tYp(new scalarField(curPatch.size()));
     scalarField& Yp = tYp();
 
-    if (typeid(curPatch) == typeid(wallFvPatch))
+    if (isType<wallFvPatch>(curPatch))
     {
-        Yp = pow(Cmu_.value(), 0.25)*y_[patchNo]
+        Yp = pow(Cmu_.value(), 0.25)
+            *y_[patchNo]
             *sqrt(k()().boundaryField()[patchNo].patchInternalField())
             /nu().boundaryField()[patchNo];
     }
@@ -158,9 +160,8 @@ tmp<scalarField> RASModel::yPlus(const label patchNo) const
     {
         WarningIn
         (
-            "tmp<scalarField> RASModel::yPlus(const label patchNo)"
-        )   << "const : " << nl
-            << "Patch " << patchNo << " is not a wall.  Returning zero field"
+            "tmp<scalarField> RASModel::yPlus(const label patchNo) const"
+        )   << "Patch " << patchNo << " is not a wall. Returning null field"
             << nl << endl;
 
         Yp.setSize(0);
@@ -185,8 +186,8 @@ bool RASModel::read()
     {
         lookup("turbulence") >> turbulence_;
         coeffDict_ = subDict(type() + "Coeffs");
-        wallFunctionDict_ = subDict("wallFunctionCoeffs");
 
+        wallFunctionDict_ = subDict("wallFunctionCoeffs");
         kappa_.readIfPresent(wallFunctionDict_);
         E_.readIfPresent(wallFunctionDict_);
         Cmu_.readIfPresent(wallFunctionDict_);

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

@@ -156,7 +156,7 @@ RNGkEpsilon::RNGkEpsilon
         autoCreateNut("nut", mesh_)
     )
 {
-    nut_ = Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    nut_ == Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
     nut_.correctBoundaryConditions();
 
     printCoeffs();
@@ -297,7 +297,7 @@ void RNGkEpsilon::correct()
 
 
     // Re-calculate viscosity
-    nut_ = Cmu_*sqr(k_)/epsilon_;
+    nut_ == Cmu_*sqr(k_)/epsilon_;
     nut_.correctBoundaryConditions();
 }
 

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

@@ -160,7 +160,7 @@ void omegaWallFunctionFvPatchScalarField::updateCoeffs()
 
         scalar yPlus = Cmu25*y[faceI]*sqrt(k[faceCellI])/nuw[faceI];
 
-        omega[faceCellI] = sqrt(k[faceCellI])/(Cmu25*kappa*y[faceCellI]);
+        omega[faceCellI] = sqrt(k[faceCellI])/(Cmu25*kappa*y[faceI]);
 
         if (yPlus > yPlusLam)
         {
@@ -168,7 +168,7 @@ void omegaWallFunctionFvPatchScalarField::updateCoeffs()
                 (nutw[faceI] + nuw[faceI])
                *magGradUw[faceI]
                *Cmu25*sqrt(k[faceCellI])
-               /(kappa*y[faceCellI]);
+               /(kappa*y[faceI]);
         }
         else
         {

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

@@ -130,7 +130,7 @@ kEpsilon::kEpsilon
         autoCreateNut("nut", mesh_)
     )
 {
-    nut_ = Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    nut_ == Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
     nut_.correctBoundaryConditions();
 
     printCoeffs();
@@ -267,7 +267,7 @@ void kEpsilon::correct()
 
 
     // Re-calculate viscosity
-    nut_ = Cmu_*sqr(k_)/epsilon_;
+    nut_ == Cmu_*sqr(k_)/epsilon_;
     nut_.correctBoundaryConditions();
 }
 

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

@@ -141,7 +141,7 @@ kOmega::kOmega
         autoCreateNut("nut", mesh_)
     )
 {
-    nut_ = k_/(omega_ + omegaSmall_);
+    nut_ == k_/(omega_ + omegaSmall_);
     nut_.correctBoundaryConditions();
 
     printCoeffs();
@@ -273,7 +273,7 @@ void kOmega::correct()
 
 
     // Re-calculate viscosity
-    nut_ = k_/omega_;
+    nut_ == k_/omega_;
     nut_.correctBoundaryConditions();
 }
 

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

@@ -250,7 +250,7 @@ kOmegaSST::kOmegaSST
         autoCreateNut("nut", mesh_)
     )
 {
-    nut_ =
+    nut_ ==
         a1_*k_
        /max
         (
@@ -411,7 +411,7 @@ void kOmegaSST::correct()
 
 
     // Re-calculate viscosity
-    nut_ = a1_*k_/max(a1_*omega_, F2()*sqrt(S2));
+    nut_ == a1_*k_/max(a1_*omega_, F2()*sqrt(S2));
     nut_.correctBoundaryConditions();
 }
 

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

@@ -182,7 +182,7 @@ realizableKE::realizableKE
     bound(k_, k0_);
     bound(epsilon_, epsilon0_);
 
-    nut_ = rCmu(fvc::grad(U_))*sqr(k_)/(epsilon_ + epsilonSmall_);
+    nut_ == rCmu(fvc::grad(U_))*sqr(k_)/(epsilon_ + epsilonSmall_);
     nut_.correctBoundaryConditions();
 
     printCoeffs();
@@ -326,7 +326,7 @@ void realizableKE::correct()
 
 
     // Re-calculate viscosity
-    nut_ = rCmu(gradU, S2, magS)*sqr(k_)/epsilon_;
+    nut_ == rCmu(gradU, S2, magS)*sqr(k_)/epsilon_;
     nut_.correctBoundaryConditions();
 }