From d1260a58e9b39b3e5015aad618859c603fc88941 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Fri, 26 Sep 2008 10:33:00 +0100
Subject: [PATCH 01/10] no pointField construction

---
 src/OpenFOAM/meshes/meshShapes/face/face.C | 25 +++++++++++++---------
 1 file changed, 15 insertions(+), 10 deletions(-)

diff --git a/src/OpenFOAM/meshes/meshShapes/face/face.C b/src/OpenFOAM/meshes/meshShapes/face/face.C
index bb0492735f..506600c31d 100644
--- a/src/OpenFOAM/meshes/meshShapes/face/face.C
+++ b/src/OpenFOAM/meshes/meshShapes/face/face.C
@@ -526,7 +526,7 @@ Foam::point Foam::face::centre(const pointField& meshPoints) const
 }
 
 
-Foam::vector Foam::face::normal(const pointField& meshPoints) const
+Foam::vector Foam::face::normal(const pointField& p) const
 {
     // Calculate the normal by summing the face triangle normals.
     // Changed to deal with small concavity by using a central decomposition
@@ -539,17 +539,22 @@ Foam::vector Foam::face::normal(const pointField& meshPoints) const
     {
         return triPointRef
         (
-            meshPoints[operator[](0)],
-            meshPoints[operator[](1)],
-            meshPoints[operator[](2)]
+            p[operator[](0)],
+            p[operator[](1)],
+            p[operator[](2)]
         ).normal();
     }
 
+    label nPoints = size();
+
     vector n = vector::zero;
 
-    point centrePoint = Foam::average(points(meshPoints));
-
-    label nPoints = size();
+    point centrePoint = vector::zero;
+    for (pI = 0; pI < nPoints; pI++)
+    {
+        centrePoint += p[operator[](pI)];
+    }
+    centrePoint /= nPoints;
 
     point nextPoint = centrePoint;
 
@@ -559,18 +564,18 @@ Foam::vector Foam::face::normal(const pointField& meshPoints) const
     {
         if (pI < nPoints - 1)
         {
-            nextPoint = meshPoints[operator[](pI + 1)];
+            nextPoint = p[operator[](pI + 1)];
         }
         else
         {
-            nextPoint = meshPoints[operator[](0)];
+            nextPoint = p[operator[](0)];
         }
 
         // Note: for best accuracy, centre point always comes last
         //
         n += triPointRef
         (
-            meshPoints[operator[](pI)],
+            p[operator[](pI)],
             nextPoint,
             centrePoint
         ).normal();

From e81c5a29e3971765e71f8ed20c99a595b0ea0299 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Fri, 26 Sep 2008 10:33:25 +0100
Subject: [PATCH 02/10] annoying message

---
 src/meshTools/sets/topoSets/faceSet.C | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/src/meshTools/sets/topoSets/faceSet.C b/src/meshTools/sets/topoSets/faceSet.C
index 7c195ca0a9..95d8dcefe9 100644
--- a/src/meshTools/sets/topoSets/faceSet.C
+++ b/src/meshTools/sets/topoSets/faceSet.C
@@ -205,11 +205,12 @@ void faceSet::sync(const polyMesh& mesh)
 
     reduce(nAdded, sumOp<label>());
 
-    if (nAdded > 0)
-    {
-        Info<< "Added an additional " << nAdded << " faces on coupled patches. "
-            << "(processorPolyPatch, cyclicPolyPatch)" << endl;
-    }
+    //if (nAdded > 0)
+    //{
+    //    Info<< "Added an additional " << nAdded
+    //        << " faces on coupled patches. "
+    //        << "(processorPolyPatch, cyclicPolyPatch)" << endl;
+    //}
 }
 
 

From cbda73951cb598a118143d1535f6b54faa094f0a Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Fri, 26 Sep 2008 10:33:46 +0100
Subject: [PATCH 03/10] regular expression wrapping

---
 src/OSspecific/Unix/regularExpression.H | 121 ++++++++++++++++++++++++
 1 file changed, 121 insertions(+)
 create mode 100644 src/OSspecific/Unix/regularExpression.H

diff --git a/src/OSspecific/Unix/regularExpression.H b/src/OSspecific/Unix/regularExpression.H
new file mode 100644
index 0000000000..dc574a4c9f
--- /dev/null
+++ b/src/OSspecific/Unix/regularExpression.H
@@ -0,0 +1,121 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::regularExpression
+
+Description
+    Wrapper around regular expressions.
+
+SourceFiles
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef regularExpression_H
+#define regularExpression_H
+
+#include <sys/types.h>
+#include <regex.h>
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class regularExpression Declaration
+\*---------------------------------------------------------------------------*/
+
+class regularExpression
+{
+    // Private data
+
+        //- Precompiled regular expression
+        regex_t* preg_;
+
+
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        regularExpression(const regularExpression&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const regularExpression&);
+
+public:
+
+
+    // Constructors
+
+        //- Construct from string
+        inline regularExpression(const string& s)
+        {
+            preg_ = new regex_t;
+
+            if (regcomp(preg_, s.c_str(), REG_EXTENDED|REG_NOSUB) != 0)
+            {
+                FatalErrorIn
+                (
+                    "regularExpression::regularExpression(const char*)"
+                )   << "Failed to compile regular expression " << s
+                    << exit(FatalError);
+            }
+        }
+
+
+    // Destructor
+
+        //- Construct from string
+        inline ~regularExpression()
+        {
+            if (preg_)
+            {
+                regfree(preg_);
+                delete preg_;
+            }
+        }
+
+
+    // Member functions
+
+        //- Matches?
+        inline bool matches(const string& s)
+        {
+            size_t nmatch = 0;
+            regmatch_t *pmatch = NULL;
+
+            return regexec(preg_, s.c_str(), nmatch, pmatch, 0) == 0;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

From 749a4bb510a3863f30fd1bbcd6d9c05d83dff20e Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Fri, 26 Sep 2008 13:10:12 +0100
Subject: [PATCH 04/10] typo

---
 src/OpenFOAM/meshes/meshShapes/face/face.C | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/OpenFOAM/meshes/meshShapes/face/face.C b/src/OpenFOAM/meshes/meshShapes/face/face.C
index 506600c31d..502e3e3eec 100644
--- a/src/OpenFOAM/meshes/meshShapes/face/face.C
+++ b/src/OpenFOAM/meshes/meshShapes/face/face.C
@@ -547,7 +547,7 @@ Foam::vector Foam::face::normal(const pointField& p) const
 
     label nPoints = size();
 
-    vector n = vector::zero;
+    register label pI;
 
     point centrePoint = vector::zero;
     for (pI = 0; pI < nPoints; pI++)
@@ -556,9 +556,9 @@ Foam::vector Foam::face::normal(const pointField& p) const
     }
     centrePoint /= nPoints;
 
-    point nextPoint = centrePoint;
+    vector n = vector::zero;
 
-    register label pI;
+    point nextPoint = centrePoint;
 
     for (pI = 0; pI < nPoints; pI++)
     {

From 1ae6914a8e519a29abb74739c6432ad8a196d161 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Fri, 26 Sep 2008 13:10:27 +0100
Subject: [PATCH 05/10] split off regular expressions

---
 .../primitives/Lists/stringListTemplates.C    | 23 +++----------------
 1 file changed, 3 insertions(+), 20 deletions(-)

diff --git a/src/OpenFOAM/primitives/Lists/stringListTemplates.C b/src/OpenFOAM/primitives/Lists/stringListTemplates.C
index f4fb22bdac..a93c46e8ac 100644
--- a/src/OpenFOAM/primitives/Lists/stringListTemplates.C
+++ b/src/OpenFOAM/primitives/Lists/stringListTemplates.C
@@ -25,9 +25,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "labelList.H"
-
-#include <sys/types.h>
-#include <regex.h>
+#include "regularExpression.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -41,24 +39,12 @@ labelList findStrings(const string& regexp, const StringList& sl)
 {
     labelList matches(sl.size());
 
-    regex_t *preg = new regex_t;
-
-    if (regcomp(preg, regexp.c_str(), REG_EXTENDED|REG_NOSUB) != 0)
-    {
-        WarningIn("findStrings(const string& regexp, const stringList& sl)")
-            << "Failed to compile regular expression " << regexp
-            << endl;
-
-        return matches;
-    }
-
-    size_t nmatch = 0;
-    regmatch_t *pmatch = NULL;
+    regularExpression re(regexp);
 
     label matchi = 0;
     forAll(sl, i)
     {
-        if (regexec(preg, sl[i].c_str(), nmatch, pmatch, 0) == 0)
+        if (re.matches(sl[i]))
         {
             matches[matchi++] = i;
         }
@@ -66,9 +52,6 @@ labelList findStrings(const string& regexp, const StringList& sl)
 
     matches.setSize(matchi);
 
-    regfree(preg);
-    delete preg;
-
     return matches;
 }
 

From e548936cb035df10b9423cef59b9b4c81cdd7551 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Fri, 26 Sep 2008 13:11:31 +0100
Subject: [PATCH 06/10] revert to pipe - behaves better

---
 wmake/wmakeScheduler | 16 +++++++++++++---
 1 file changed, 13 insertions(+), 3 deletions(-)

diff --git a/wmake/wmakeScheduler b/wmake/wmakeScheduler
index bf9bfc914b..2a2ed0fd76 100755
--- a/wmake/wmakeScheduler
+++ b/wmake/wmakeScheduler
@@ -137,6 +137,16 @@ done
 # is returned and not of colouring pipe.
 set -o pipefail
 
+# Define function to colour output by argument 1
+colourPipe(){
+    if [ "$1" ]; then
+        (while read line; do setterm -foreground $1; echo "$line" ; done; setterm -foreground default)
+    else
+        cat
+    fi
+}
+
+
 colourIndex=0
 
 while :
@@ -156,14 +166,14 @@ do
             if lockfile -r0 "$lockFile" 2>/dev/null; then
                 if [ "$WM_COLOURS" ]; then
                     # Set colour
-                    colourString=`setterm -foreground ${colours[$colourIndex]}`
+                    colour="${colours[$colourIndex]}"
 
                     if [ "$host" = "$HOST" ]; then
-                        eval $* 2>&1 | sed -e "s/^/$colourString/"
+                        eval $* 2>&1 | colourPipe "$colour"
                     elif [ -n "$JOB_ID" ]; then
                         qrsh -inherit -v PWD $host "$rcmd"
                     else
-                        ssh $host "$sourceFoam 2>/dev/null; cd $PWD && $rcmd" 2>&1 | sed -e "s/^/$colourString/"
+                        ssh $host "$sourceFoam 2>/dev/null; cd $PWD && $rcmd" 2>&1 | colourPipe "$colour"
                     fi
                     retval=$?
                 else

From 1bf595225dc2a9cb2d5df6d5df9a3f493eae66a1 Mon Sep 17 00:00:00 2001
From: andy <a.heather@opencfd.co.uk>
Date: Fri, 26 Sep 2008 14:29:17 +0100
Subject: [PATCH 07/10] adding minMaxFields

---
 src/postProcessing/Allwmake | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/postProcessing/Allwmake b/src/postProcessing/Allwmake
index f67ba57aa1..1d889e8b76 100755
--- a/src/postProcessing/Allwmake
+++ b/src/postProcessing/Allwmake
@@ -6,5 +6,6 @@ wmake libo postCalc
 wmake libso forces
 wmake libso fieldAverage
 wmake libso foamCalcFunctions
+wmake libso minMaxFields
 
 # ----------------------------------------------------------------- end-of-file

From 26f2c24d84d653169811d752b5ff3884c2c9c87a Mon Sep 17 00:00:00 2001
From: henry <Henry Weller h.weller@opencfd.co.uk>
Date: Sat, 27 Sep 2008 09:25:30 +0100
Subject: [PATCH 08/10] Added support for quadraticFit interpolation scheme.
 Added the quadraticFit interpolation scheme.

---
 src/OpenFOAM/Make/files                       |   1 +
 .../matrices/DiagonalMatrix/DiagonalMatrix.C  |  99 ++++
 .../matrices/DiagonalMatrix/DiagonalMatrix.H  | 102 ++++
 .../LUscalarMatrix/LUscalarMatrixTemplates.C  |  14 +-
 .../{containers => matrices}/Matrix/Matrix.C  | 174 +++---
 .../{containers => matrices}/Matrix/Matrix.H  |  83 ++-
 .../{containers => matrices}/Matrix/MatrixI.H |  52 +-
 .../Matrix/MatrixIO.C                         |  44 +-
 src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.C  | 400 +++++++++++++
 src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.H  | 128 +++++
 src/OpenFOAM/matrices/scalarMatrix/SVD/SVDI.H |  75 +++
 .../matrices/scalarMatrix/scalarMatrix.C      | 159 +++++-
 .../matrices/scalarMatrix/scalarMatrix.H      |  46 +-
 .../scalarMatrix/scalarMatrixTemplates.C      |  26 +-
 .../matrices/simpleMatrix/simpleMatrix.C      |  68 +--
 .../matrices/simpleMatrix/simpleMatrix.H      |  66 +--
 src/finiteVolume/Make/files                   |   5 +
 .../fvMesh/extendedStencil/extendedStencil.C  | 525 ++++++++++++++++++
 .../fvMesh/extendedStencil/extendedStencil.H  | 151 +++++
 .../extendedStencilTemplates.C                | 129 +++++
 .../schemes/quadraticFit/quadraticFit.C       |  36 ++
 .../schemes/quadraticFit/quadraticFit.H       | 138 +++++
 .../schemes/quadraticFit/quadraticFitData.C   | 310 +++++++++++
 .../schemes/quadraticFit/quadraticFitData.H   | 140 +++++
 24 files changed, 2722 insertions(+), 249 deletions(-)
 create mode 100644 src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C
 create mode 100644 src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H
 rename src/OpenFOAM/{containers => matrices}/Matrix/Matrix.C (62%)
 rename src/OpenFOAM/{containers => matrices}/Matrix/Matrix.H (70%)
 rename src/OpenFOAM/{containers => matrices}/Matrix/MatrixI.H (66%)
 rename src/OpenFOAM/{containers => matrices}/Matrix/MatrixIO.C (84%)
 create mode 100644 src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.C
 create mode 100644 src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.H
 create mode 100644 src/OpenFOAM/matrices/scalarMatrix/SVD/SVDI.H
 create mode 100644 src/finiteVolume/fvMesh/extendedStencil/extendedStencil.C
 create mode 100644 src/finiteVolume/fvMesh/extendedStencil/extendedStencil.H
 create mode 100644 src/finiteVolume/fvMesh/extendedStencil/extendedStencilTemplates.C
 create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.C
 create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.H
 create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C
 create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.H

diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files
index 7b137dedc0..57dc7f8dc7 100644
--- a/src/OpenFOAM/Make/files
+++ b/src/OpenFOAM/Make/files
@@ -179,6 +179,7 @@ matrices/solution/solution.C
 
 scalarMatrix = matrices/scalarMatrix
 $(scalarMatrix)/scalarMatrix.C
+$(scalarMatrix)/SVD/SVD.C
 
 LUscalarMatrix = matrices/LUscalarMatrix
 $(LUscalarMatrix)/LUscalarMatrix.C
diff --git a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C
new file mode 100644
index 0000000000..c023e59c5a
--- /dev/null
+++ b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C
@@ -0,0 +1,99 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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>
+Foam::DiagonalMatrix<Type>::DiagonalMatrix(const Matrix<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;
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H
new file mode 100644
index 0000000000..36138fc25d
--- /dev/null
+++ b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H
@@ -0,0 +1,102 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 Type> class Matrix;
+
+/*---------------------------------------------------------------------------*\
+                           Class DiagonalMatrix Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class DiagonalMatrix
+:
+    public List<Type>
+{
+public:
+
+    // Constructors
+
+        //- Construct from diagonal component of a Matrix
+        DiagonalMatrix<Type>(const Matrix<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
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C
index cbb59e2c86..ac0bc17207 100644
--- a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C
+++ b/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)
                 );
             }
         }
diff --git a/src/OpenFOAM/containers/Matrix/Matrix.C b/src/OpenFOAM/matrices/Matrix/Matrix.C
similarity index 62%
rename from src/OpenFOAM/containers/Matrix/Matrix.C
rename to src/OpenFOAM/matrices/Matrix/Matrix.C
index 733e04034e..4fbef71c2b 100644
--- a/src/OpenFOAM/containers/Matrix/Matrix.C
+++ b/src/OpenFOAM/matrices/Matrix/Matrix.C
@@ -28,13 +28,13 @@ License
 
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
-template<class T>
-void Foam::Matrix<T>::allocate()
+template<class Type>
+void Foam::Matrix<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 Type>
+Foam::Matrix<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 Type>
+const Foam::Matrix<Type>& Foam::Matrix<Type>::null()
 {
-    Matrix<T>* nullPtr = reinterpret_cast<Matrix<T>*>(NULL);
+    Matrix<Type>* nullPtr = reinterpret_cast<Matrix<Type>*>(NULL);
     return *nullPtr;
 }
 
 
-template<class T>
-Foam::Matrix<T>::Matrix(const label n, const label m)
+template<class Type>
+Foam::Matrix<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<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 Type>
+Foam::Matrix<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<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 Type>
+Foam::Matrix<Type>::Matrix(const Matrix<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 Type>
+void Foam::Matrix<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 Type>
+void Foam::Matrix<Type>::transfer(Matrix<Type>& a)
 {
     clear();
 
@@ -171,12 +171,30 @@ void Foam::Matrix<T>::transfer(Matrix<T>& a)
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
-template<class T>
-void Foam::Matrix<T>::operator=(const T& t)
+template<class Type>
+Foam::Matrix<Type> Foam::Matrix<Type>::T() const
+{
+    const Matrix<Type>& A = *this;
+    Matrix<Type> 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;
+}
+
+
+template<class Type>
+void Foam::Matrix<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 Type>
+void Foam::Matrix<Type>::operator=(const Matrix<Type>& a)
 {
     if (this == &a)
     {
-        FatalErrorIn("Matrix<T>::operator=(const Matrix<T>&)")
+        FatalErrorIn("Matrix<Type>::operator=(const Matrix<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 Type>
+const Type& Foam::max(const Matrix<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<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 Type>
+const Type& Foam::min(const Matrix<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<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 Type>
+Foam::Matrix<Type> Foam::operator-(const Matrix<Type>& a)
 {
-    Matrix<T> na(a.n_, a.m_);
+    Matrix<Type> 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 Type>
+Foam::Matrix<Type> Foam::operator+(const Matrix<Type>& a, const Matrix<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<Type>::operator+(const Matrix<Type>&, const Matrix<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<Type>::operator+(const Matrix<Type>&, const Matrix<Type>&)"
+        )   << "attempted add matrices with different number of columns: "
             << a.m_ << ", " << b.m_
             << abort(FatalError);
     }
 
-    Matrix<T> ab(a.n_, a.m_);
+    Matrix<Type> 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 Type>
+Foam::Matrix<Type> Foam::operator-(const Matrix<Type>& a, const Matrix<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<Type>::operator-(const Matrix<Type>&, const Matrix<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<Type>::operator-(const Matrix<Type>&, const Matrix<Type>&)"
+        )   << "attempted add matrices with different number of columns: "
             << a.m_ << ", " << b.m_
             << abort(FatalError);
     }
 
-    Matrix<T> ab(a.n_, a.m_);
+    Matrix<Type> 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,15 +405,15 @@ 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 Type>
+Foam::Matrix<Type> Foam::operator*(const scalar s, const Matrix<Type>& a)
 {
-    Matrix<T> sa(a.n_, a.m_);
+    Matrix<Type> 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_;;
         for (register label i=0; i<nm; i++)
diff --git a/src/OpenFOAM/containers/Matrix/Matrix.H b/src/OpenFOAM/matrices/Matrix/Matrix.H
similarity index 70%
rename from src/OpenFOAM/containers/Matrix/Matrix.H
rename to src/OpenFOAM/matrices/Matrix/Matrix.H
index ab68ce5bad..6b4da39d4f 100644
--- a/src/OpenFOAM/containers/Matrix/Matrix.H
+++ b/src/OpenFOAM/matrices/Matrix/Matrix.H
@@ -51,25 +51,37 @@ namespace Foam
 
 // Forward declaration of friend functions and operators
 
-template<class T> class Matrix;
+template<class Type> class Matrix;
 
-template<class T> const T& max(const Matrix<T>&);
-template<class T> const T& min(const Matrix<T>&);
+template<class Type> const Type& max(const Matrix<Type>&);
+template<class Type> const Type& min(const Matrix<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 Type> Matrix<Type> operator-(const Matrix<Type>&);
+template<class Type> Matrix<Type> operator+
+(
+    const Matrix<Type>&,
+    const Matrix<Type>&
+);
+template<class Type> Matrix<Type> operator-
+(
+    const Matrix<Type>&,
+    const Matrix<Type>&
+);
+template<class Type> Matrix<Type> operator*
+(
+    const scalar,
+    const Matrix<Type>&
+);
 
-template<class T> Istream& operator>>(Istream&, Matrix<T>&);
-template<class T> Ostream& operator<<(Ostream&, const Matrix<T>&);
+template<class Type> Istream& operator>>(Istream&, Matrix<Type>&);
+template<class Type> Ostream& operator<<(Ostream&, const Matrix<Type>&);
 
 
 /*---------------------------------------------------------------------------*\
                            Class Matrix Declaration
 \*---------------------------------------------------------------------------*/
 
-template<class T>
+template<class Type>
 class Matrix
 {
     // Private data
@@ -78,7 +90,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 +109,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<Type>&);
 
         //- Construct from Istream.
         Matrix(Istream&);
 
         //- Clone
-        inline autoPtr<Matrix<T> > clone() const;
+        inline autoPtr<Matrix<Type> > clone() const;
 
 
     // Destructor
@@ -117,7 +129,7 @@ public:
     // Member functions
 
         //- Return a null Matrix
-        static const Matrix<T>& null();
+        static const Matrix<Type>& null();
 
 
         // Access
@@ -148,45 +160,60 @@ public:
 
             //- Transfer the contents of the argument Matrix into this Matrix
             //  and annull the argument Matrix.
-            void transfer(Matrix<T>&);
+            void transfer(Matrix<Type>&);
 
 
     // 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;
+
+        //- Return the transpose of the matrix
+        Matrix<Type> T() const;
 
         //- Assignment operator. Takes linear time.
-        void operator=(const Matrix<T>&);
+        void operator=(const Matrix<Type>&);
 
         //- Assignment of all entries to the given value
-        void operator=(const T&);
+        void operator=(const Type&);
 
 
     // Friend functions
 
-        friend const T& max<T>(const Matrix<T>&);
-        friend const T& min<T>(const Matrix<T>&);
+        friend const Type& max<Type>(const Matrix<Type>&);
+        friend const Type& min<Type>(const Matrix<Type>&);
 
 
     // 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>&);
+        friend Matrix<Type> operator-<Type>(const Matrix<Type>&);
+        friend Matrix<Type> operator+<Type>
+        (
+            const Matrix<Type>&,
+            const Matrix<Type>&
+        );
+        friend Matrix<Type> operator-<Type>
+        (
+            const Matrix<Type>&,
+            const Matrix<Type>&
+        );
+        friend Matrix<Type> operator*<Type>
+        (
+            const scalar,
+            const Matrix<Type>&
+        );
 
 
     // IOstream operators
 
         //- Read Matrix from Istream, discarding contents of existing Matrix.
-        friend Istream& operator>> <T>(Istream&, Matrix<T>&);
+        friend Istream& operator>> <Type>(Istream&, Matrix<Type>&);
 
         // Write Matrix to Ostream.
-        friend Ostream& operator<< <T>(Ostream&, const Matrix<T>&);
+        friend Ostream& operator<< <Type>(Ostream&, const Matrix<Type>&);
 };
 
 
diff --git a/src/OpenFOAM/containers/Matrix/MatrixI.H b/src/OpenFOAM/matrices/Matrix/MatrixI.H
similarity index 66%
rename from src/OpenFOAM/containers/Matrix/MatrixI.H
rename to src/OpenFOAM/matrices/Matrix/MatrixI.H
index d0e3fe8fdf..8323174d4f 100644
--- a/src/OpenFOAM/containers/Matrix/MatrixI.H
+++ b/src/OpenFOAM/matrices/Matrix/MatrixI.H
@@ -26,8 +26,8 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-template<class T>
-inline Foam::Matrix<T>::Matrix()
+template<class Type>
+inline Foam::Matrix<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 Type>
+inline Foam::autoPtr<Foam::Matrix<Type> > Foam::Matrix<Type>::clone() const
 {
-    return autoPtr<Matrix<T> >(new Matrix<T>(*this));
+    return autoPtr<Matrix<Type> >(new Matrix<Type>(*this));
 }
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 //- Return the number of rows
-template<class T>
-inline Foam::label Foam::Matrix<T>::n() const
+template<class Type>
+inline Foam::label Foam::Matrix<Type>::n() const
 {
     return n_;
 }
 
 
-//- Return the number of columns
-template<class T>
-inline Foam::label Foam::Matrix<T>::m() const
+template<class Type>
+inline Foam::label Foam::Matrix<Type>::m() const
 {
     return m_;
 }
 
 
-//- Return the number of columns
-template<class T>
-inline Foam::label Foam::Matrix<T>::size() const
+template<class Type>
+inline Foam::label Foam::Matrix<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 Type>
+inline void Foam::Matrix<Type>::checki(const label i) const
 {
     if (!n_)
     {
-        FatalErrorIn("Matrix<T>::checki(const label)")
+        FatalErrorIn("Matrix<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<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 Type>
+inline void Foam::Matrix<Type>::checkj(const label j) const
 {
     if (!m_)
     {
-        FatalErrorIn("Matrix<T>::checkj(const label)")
+        FatalErrorIn("Matrix<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<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 Type>
+inline Type* Foam::Matrix<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 Type>
+inline const Type* Foam::Matrix<Type>::operator[](const label i) const
 {
 #   ifdef FULLDEBUG
     checki(i);
diff --git a/src/OpenFOAM/containers/Matrix/MatrixIO.C b/src/OpenFOAM/matrices/Matrix/MatrixIO.C
similarity index 84%
rename from src/OpenFOAM/containers/Matrix/MatrixIO.C
rename to src/OpenFOAM/matrices/Matrix/MatrixIO.C
index 710e42876e..9b5a922f0d 100644
--- a/src/OpenFOAM/containers/Matrix/MatrixIO.C
+++ b/src/OpenFOAM/matrices/Matrix/MatrixIO.C
@@ -32,8 +32,8 @@ License
 
 // * * * * * * * * * * * * * * * Ostream Operator *  * * * * * * * * * * * * //
 
-template<class T>
-Foam::Matrix<T>::Matrix(Istream& is)
+template<class Type>
+Foam::Matrix<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 Type>
+Foam::Istream& Foam::operator>>(Istream& is, Matrix<Type>& M)
 {
     // Anull matrix
     M.clear();
 
-    is.fatalCheck("operator>>(Istream&, Matrix<T>&)");
+    is.fatalCheck("operator>>(Istream&, Matrix<Type>&)");
 
     token firstToken(is);
 
-    is.fatalCheck("operator>>(Istream&, Matrix<T>&) : reading first token");
+    is.fatalCheck("operator>>(Istream&, Matrix<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<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<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<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<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 Type>
+Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<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));
         }
     }
 
diff --git a/src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.C b/src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.C
new file mode 100644
index 0000000000..a1355c393d
--- /dev/null
+++ b/src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.C
@@ -0,0 +1,400 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "scalarMatrix.H"
+#include "ListOps.H"
+
+// * * * * * * * * * * * * * * * * Constructor  * * * * * * * * * * * * * * //
+
+Foam::SVD::SVD(const Matrix<scalar>& 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(Matrix<scalar>& A)")
+                    << "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
+    /*Matrix<scalar> 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;
+    }
+    */
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.H b/src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.H
new file mode 100644
index 0000000000..b597027732
--- /dev/null
+++ b/src/OpenFOAM/matrices/scalarMatrix/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
+    SVD
+
+Description
+    Singular value decomposition of a rectangular matrix.
+
+SourceFiles
+    SVD.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef SVD_H
+#define SVD_H
+
+#include "DiagonalMatrix.H"
+#include "Matrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+/*---------------------------------------------------------------------------*\
+                      Class SVD Declaration
+\*---------------------------------------------------------------------------*/
+
+class SVD
+{
+    // Private data
+
+        //- Rectangular matrix with the same dimensions as the input
+        Matrix<scalar> U_;
+
+        //- square matrix V
+        Matrix<scalar> V_;
+
+        //- The singular values
+        DiagonalMatrix<scalar> S_;
+
+        //- The matrix product V S^(-1) U^T
+        Matrix<scalar> 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 Matrix<scalar>& A, const scalar minCondition = 0);
+
+
+    // Access functions
+
+        //- Return U
+        inline const Matrix<scalar>& U() const;
+
+        //- Return the square matrix V
+        inline const Matrix<scalar>& V() const;
+
+        //- Return the singular values
+        inline const DiagonalMatrix<scalar>& S() const;
+
+        //- Return VSinvUt (the pseudo inverse)
+        inline const Matrix<scalar>& 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
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/scalarMatrix/SVD/SVDI.H b/src/OpenFOAM/matrices/scalarMatrix/SVD/SVDI.H
new file mode 100644
index 0000000000..337647edfb
--- /dev/null
+++ b/src/OpenFOAM/matrices/scalarMatrix/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::Matrix<Foam::scalar>& Foam::SVD::U() const
+{
+    return U_;
+}
+
+inline const Foam::Matrix<Foam::scalar>& Foam::SVD::V() const
+{
+    return V_;
+}
+
+inline const Foam::DiagonalMatrix<Foam::scalar>& Foam::SVD::S() const
+{
+    return S_;
+}
+
+inline const Foam::Matrix<Foam::scalar>& 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;
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.C b/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.C
index 2ff9e39f18..b88e167b6f 100644
--- a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.C
+++ b/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.C
@@ -25,6 +25,8 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "scalarMatrix.H"
+#include "SVD.H"
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -136,7 +138,7 @@ void Foam::scalarMatrix::LUDecompose
             {
                 Swap(matrixj[k], matrixiMax[k]);
             }
-            
+
             vv[iMax] = vv[j];
         }
 
@@ -158,4 +160,159 @@ void Foam::scalarMatrix::LUDecompose
 }
 
 
+// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
+
+void Foam::multiply
+(
+    Matrix<scalar>& ans,         // value changed in return
+    const Matrix<scalar>& A,
+    const Matrix<scalar>& B
+)
+{
+    if (A.m() != B.n())
+    {
+        FatalErrorIn
+        (
+            "multiply("
+            "Matrix<scalar>& answer "
+            "const Matrix<scalar>& A, "
+            "const Matrix<scalar>& B)"
+        )   << "A and B must have identical inner dimensions but A.m = "
+            << A.m() << " and B.n = " << B.n()
+            << abort(FatalError);
+    }
+
+    ans = Matrix<scalar>(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
+(
+    Matrix<scalar>& ans,         // value changed in return
+    const Matrix<scalar>& A,
+    const Matrix<scalar>& B,
+    const Matrix<scalar>& C
+)
+{
+    if (A.m() != B.n())
+    {
+        FatalErrorIn
+        (
+            "multiply("
+            "const Matrix<scalar>& A, "
+            "const Matrix<scalar>& B, "
+            "const Matrix<scalar>& C, "
+            "Matrix<scalar>& 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 Matrix<scalar>& A, "
+            "const Matrix<scalar>& B, "
+            "const Matrix<scalar>& C, "
+            "Matrix<scalar>& answer)"
+        )   << "B and C must have identical inner dimensions but B.m = "
+            << B.m() << " and C.n = " << C.n()
+            << abort(FatalError);
+    }
+
+    ans = Matrix<scalar>(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
+(
+    Matrix<scalar>& ans,         // value changed in return
+    const Matrix<scalar>& A,
+    const DiagonalMatrix<scalar>& B,
+    const Matrix<scalar>& C
+)
+{
+    if (A.m() != B.size())
+    {
+        FatalErrorIn
+        (
+            "multiply("
+            "const Matrix<scalar>& A, "
+            "const DiagonalMatrix<scalar>& B, "
+            "const Matrix<scalar>& C, "
+            "Matrix<scalar>& 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 Matrix<scalar>& A, "
+            "const DiagonalMatrix<scalar>& B, "
+            "const Matrix<scalar>& C, "
+            "Matrix<scalar>& answer)"
+        )   << "B and C must have identical inner dimensions but B.m = "
+            << B.size() << " and C.n = " << C.n()
+            << abort(FatalError);
+    }
+
+    ans = Matrix<scalar>(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::Matrix<Foam::scalar> Foam::SVDinv
+(
+    const Matrix<scalar>& A,
+    scalar minCondition
+)
+{
+    SVD svd(A, minCondition);
+    return svd.VSinvUt();
+}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.H b/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.H
index f2b0764f1a..2d22e71d4e 100644
--- a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.H
+++ b/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.H
@@ -37,6 +37,7 @@ SourceFiles
 #define scalarMatrix_H
 
 #include "Matrix.H"
+#include "DiagonalMatrix.H"
 #include "scalarField.H"
 #include "labelList.H"
 
@@ -75,13 +76,13 @@ public:
 
         //- 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);
+        template<class Type>
+        static void solve(Matrix<scalar>& matrix, Field<Type>& 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;
+        template<class Type>
+        void solve(Field<Type>& psi, const Field<Type>& source) const;
 
 
         //- LU decompose the matrix with pivoting
@@ -93,21 +94,50 @@ public:
 
         //- LU back-substitution with given source, returning the solution
         //  in the source
-        template<class T>
+        template<class Type>
         static void LUBacksubstitute
         (
             const Matrix<scalar>& luMmatrix,
             const labelList& pivotIndices,
-            Field<T>& source
+            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 T>
-        static void LUsolve(Matrix<scalar>& matrix, Field<T>& source);
+        template<class Type>
+        static void LUsolve(Matrix<scalar>& matrix, Field<Type>& source);
 };
 
 
+// Global functions
+
+void multiply
+(
+    Matrix<scalar>& answer,         // value changed in return
+    const Matrix<scalar>& A,
+    const Matrix<scalar>& B
+);
+
+void multiply
+(
+    Matrix<scalar>& answer,         // value changed in return
+    const Matrix<scalar>& A,
+    const Matrix<scalar>& B,
+    const Matrix<scalar>& C
+);
+
+void multiply
+(
+    Matrix<scalar>& answer,         // value changed in return
+    const Matrix<scalar>& A,
+    const DiagonalMatrix<scalar>& B,
+    const Matrix<scalar>& C
+);
+
+//- Return the inverse of matrix A using SVD
+Matrix<scalar> SVDinv(const Matrix<scalar>& A, scalar minCondition = 0);
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrixTemplates.C b/src/OpenFOAM/matrices/scalarMatrix/scalarMatrixTemplates.C
index b729e32fcc..fb48a9b26d 100644
--- a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrixTemplates.C
+++ b/src/OpenFOAM/matrices/scalarMatrix/scalarMatrixTemplates.C
@@ -29,11 +29,11 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class T>
+template<class Type>
 void Foam::scalarMatrix::solve
 (
     Matrix<scalar>& tmpMatrix,
-    Field<T>& sourceSol
+    Field<Type>& sourceSol
 )
 {
     label n = tmpMatrix.n();
@@ -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,8 +101,8 @@ void Foam::scalarMatrix::solve
 }
 
 
-template<class T>
-void Foam::scalarMatrix::solve(Field<T>& psi, const Field<T>& source) const
+template<class Type>
+void Foam::scalarMatrix::solve(Field<Type>& psi, const Field<Type>& source) const
 {
     Matrix<scalar> tmpMatrix = *this;
     psi = source;
@@ -110,12 +110,12 @@ void Foam::scalarMatrix::solve(Field<T>& psi, const Field<T>& source) const
 }
 
 
-template<class T>
+template<class Type>
 void Foam::scalarMatrix::LUBacksubstitute
 (
     const Matrix<scalar>& luMatrix,
     const labelList& pivotIndices,
-    Field<T>& sourceSol
+    Field<Type>& sourceSol
 )
 {
     label n = luMatrix.n();
@@ -125,7 +125,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 +136,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 +146,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 +159,11 @@ void Foam::scalarMatrix::LUBacksubstitute
 }
 
 
-template<class T>
+template<class Type>
 void Foam::scalarMatrix::LUsolve
 (
     Matrix<scalar>& matrix,
-    Field<T>& sourceSol
+    Field<Type>& sourceSol
 )
 {
     labelList pivotIndices(matrix.n());
diff --git a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C
index 22de758f3a..1265b87ee5 100644
--- a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C
+++ b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C
@@ -28,19 +28,19 @@ 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)
+    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 Field<Type>& source
 )
 :
     scalarMatrix(matrix),
@@ -48,8 +48,8 @@ Foam::simpleMatrix<T>::simpleMatrix
 {}
 
 
-template<class T>
-Foam::simpleMatrix<T>::simpleMatrix(Istream& is)
+template<class Type>
+Foam::simpleMatrix<Type>::simpleMatrix(Istream& is)
 :
     scalarMatrix(is),
     source_(is)
@@ -58,11 +58,11 @@ Foam::simpleMatrix<T>::simpleMatrix(Istream& 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_;
+    Field<Type> sourceSol = source_;
 
     scalarMatrix::solve(tmpMatrix, sourceSol);
 
@@ -70,11 +70,11 @@ Foam::Field<T> Foam::simpleMatrix<T>::solve() const
 }
 
 
-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_;
+    Field<Type> sourceSol = source_;
 
     scalarMatrix::LUsolve(luMatrix, sourceSol);
 
@@ -84,26 +84,26 @@ 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);
     }
@@ -115,14 +115,14 @@ void Foam::simpleMatrix<T>::operator=(const simpleMatrix<T>& m)
 
 // * * * * * * * * * * * * * * * 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),
@@ -131,14 +131,14 @@ Foam::simpleMatrix<T> Foam::operator+
 }
 
 
-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),
@@ -147,17 +147,17 @@ Foam::simpleMatrix<T> Foam::operator-
 }
 
 
-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_;
     return os;
diff --git a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H
index d1c411452b..d7b1a027be 100644
--- a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H
+++ b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H
@@ -45,35 +45,35 @@ namespace Foam
 
 // Forward declaration of friend functions and operators
 
-template<class T>
+template<class Type>
 class simpleMatrix;
 
-template<class T>
-simpleMatrix<T> operator+
+template<class Type>
+simpleMatrix<Type> operator+
 (
-    const simpleMatrix<T>&,
-    const simpleMatrix<T>&
+    const simpleMatrix<Type>&,
+    const simpleMatrix<Type>&
 );
 
-template<class T>
-simpleMatrix<T> operator-
+template<class Type>
+simpleMatrix<Type> operator-
 (
-    const simpleMatrix<T>&,
-    const simpleMatrix<T>&
+    const simpleMatrix<Type>&,
+    const simpleMatrix<Type>&
 );
 
-template<class T>
-simpleMatrix<T> operator*
+template<class Type>
+simpleMatrix<Type> operator*
 (
     const scalar,
-    const simpleMatrix<T>&
+    const simpleMatrix<Type>&
 );
 
-template<class T>
+template<class Type>
 Ostream& operator<<
 (
     Ostream&,
-    const simpleMatrix<T>&
+    const simpleMatrix<Type>&
 );
 
 
@@ -81,14 +81,14 @@ Ostream& operator<<
                            Class simpleMatrix Declaration
 \*---------------------------------------------------------------------------*/
 
-template<class T>
+template<class Type>
 class simpleMatrix
 :
     public scalarMatrix
 {
     // Private data
 
-        Field<T> source_;
+        Field<Type> source_;
 
 
 public:
@@ -99,25 +99,25 @@ public:
         simpleMatrix(const label);
 
         //- Construct from components
-        simpleMatrix(const scalarMatrix&, const Field<T>&);
+        simpleMatrix(const scalarMatrix&, 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,45 +125,45 @@ 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>&);
+        void operator=(const simpleMatrix<Type>&);
 
 
     // Friend Operators
 
-        friend simpleMatrix<T> operator+ <T>
+        friend simpleMatrix<Type> operator+ <Type>
         (
-            const simpleMatrix<T>&,
-            const simpleMatrix<T>&
+            const simpleMatrix<Type>&,
+            const simpleMatrix<Type>&
         );
 
-        friend simpleMatrix<T> operator- <T>
+        friend simpleMatrix<Type> operator- <Type>
         (
-            const simpleMatrix<T>&,
-            const simpleMatrix<T>&
+            const simpleMatrix<Type>&,
+            const simpleMatrix<Type>&
         );
 
-        friend simpleMatrix<T> operator* <T>
+        friend simpleMatrix<Type> operator* <Type>
         (
             const scalar,
-            const simpleMatrix<T>&
+            const simpleMatrix<Type>&
         );
 
 
     // Ostream Operator
 
-        friend Ostream& operator<< <T>
+        friend Ostream& operator<< <Type>
         (
             Ostream&,
-            const simpleMatrix<T>&
+            const simpleMatrix<Type>&
         );
 };
 
diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index 0c3fae1887..25137c58c7 100644
--- a/src/finiteVolume/Make/files
+++ b/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
diff --git a/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.C b/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.C
new file mode 100644
index 0000000000..bb4c6c1534
--- /dev/null
+++ b/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);
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.H b/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.H
new file mode 100644
index 0000000000..1d00972403
--- /dev/null
+++ b/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
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/fvMesh/extendedStencil/extendedStencilTemplates.C b/src/finiteVolume/fvMesh/extendedStencil/extendedStencilTemplates.C
new file mode 100644
index 0000000000..a408621dc4
--- /dev/null
+++ b/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;
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.C
new file mode 100644
index 0000000000..a215576330
--- /dev/null
+++ b/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);
+}
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.H
new file mode 100644
index 0000000000..412b454015
--- /dev/null
+++ b/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
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C
new file mode 100644
index 0000000000..fd0e55fcdd
--- /dev/null
+++ b/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
+    Matrix<scalar> 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;
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.H
new file mode 100644
index 0000000000..19859ea621
--- /dev/null
+++ b/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
+
+// ************************************************************************* //

From b722041fff3129f208eccde20cac205556be1002 Mon Sep 17 00:00:00 2001
From: henry <Henry Weller h.weller@opencfd.co.uk>
Date: Sun, 28 Sep 2008 22:46:33 +0100
Subject: [PATCH 09/10] Corrected reading of "LESProperties".

---
 .../execFlowFunctionObjects/execFlowFunctionObjects.C           | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/applications/utilities/postProcessing/miscellaneous/execFlowFunctionObjects/execFlowFunctionObjects.C b/applications/utilities/postProcessing/miscellaneous/execFlowFunctionObjects/execFlowFunctionObjects.C
index b848e330ea..3f7a99040b 100644
--- a/applications/utilities/postProcessing/miscellaneous/execFlowFunctionObjects/execFlowFunctionObjects.C
+++ b/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,

From 318b71e206e9e635aa649938e9d1c5dfbc0b7f6b Mon Sep 17 00:00:00 2001
From: henry <Henry Weller h.weller@opencfd.co.uk>
Date: Sun, 28 Sep 2008 22:50:57 +0100
Subject: [PATCH 10/10] Made Matrix generic base-class for SquareMatrix and
 RectangularMatrix by introducing the type it will become as a template
 argument.

Brought everything into line with this change.
---
 src/ODE/ODE/ODE.H                             |   4 +-
 src/ODE/ODESolvers/KRR4/KRR4.C                |  19 +--
 src/ODE/ODESolvers/KRR4/KRR4.H                |   6 +-
 src/ODE/ODESolvers/SIBS/SIBS.H                |  10 +-
 src/ODE/ODESolvers/SIBS/SIMPR.C               |  13 +-
 src/ODE/ODESolvers/SIBS/polyExtrapolate.C     |   2 +-
 src/OpenFOAM/Make/files                       |   6 +-
 .../IOstreams/Pstreams/PstreamCommsStruct.C   |   2 -
 .../matrices/DiagonalMatrix/DiagonalMatrix.C  |   3 +-
 .../matrices/DiagonalMatrix/DiagonalMatrix.H  |   5 +-
 .../matrices/LUscalarMatrix/LUscalarMatrix.C  |   8 +-
 .../matrices/LUscalarMatrix/LUscalarMatrix.H  |  10 +-
 src/OpenFOAM/matrices/Matrix/Matrix.C         | 106 ++++++------
 src/OpenFOAM/matrices/Matrix/Matrix.H         | 104 ++++++------
 src/OpenFOAM/matrices/Matrix/MatrixI.H        |  46 +++---
 src/OpenFOAM/matrices/Matrix/MatrixIO.C       |  24 +--
 .../RectangularMatrix/RectangularMatrix.H     |  92 +++++++++++
 .../RectangularMatrix/RectangularMatrixI.H    |  70 ++++++++
 .../matrices/SquareMatrix/SquareMatrix.H      |  96 +++++++++++
 .../matrices/SquareMatrix/SquareMatrixI.H     |  75 +++++++++
 .../SVD/SVD.C                                 |  13 +-
 .../SVD/SVD.H                                 |  22 +--
 .../SVD/SVDI.H                                |   8 +-
 .../scalarMatrices.C}                         |  99 +++++------
 .../matrices/scalarMatrices/scalarMatrices.H  | 137 ++++++++++++++++
 .../scalarMatricesTemplates.C}                |  25 +--
 .../matrices/scalarMatrix/scalarMatrix.H      | 155 ------------------
 .../matrices/simpleMatrix/simpleMatrix.C      |  28 ++--
 .../matrices/simpleMatrix/simpleMatrix.H      |  72 +++-----
 .../schemes/quadraticFit/quadraticFitData.C   |  14 +-
 .../chemistryModel/chemistryModel.C           |  38 ++---
 .../chemistryModel/chemistryModel.H           |   3 +-
 32 files changed, 791 insertions(+), 524 deletions(-)
 create mode 100644 src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H
 create mode 100644 src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H
 create mode 100644 src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H
 create mode 100644 src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H
 rename src/OpenFOAM/matrices/{scalarMatrix => scalarMatrices}/SVD/SVD.C (96%)
 rename src/OpenFOAM/matrices/{scalarMatrix => scalarMatrices}/SVD/SVD.H (86%)
 rename src/OpenFOAM/matrices/{scalarMatrix => scalarMatrices}/SVD/SVDI.H (88%)
 rename src/OpenFOAM/matrices/{scalarMatrix/scalarMatrix.C => scalarMatrices/scalarMatrices.C} (76%)
 create mode 100644 src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H
 rename src/OpenFOAM/matrices/{scalarMatrix/scalarMatrixTemplates.C => scalarMatrices/scalarMatricesTemplates.C} (90%)
 delete mode 100644 src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.H

diff --git a/src/ODE/ODE/ODE.H b/src/ODE/ODE/ODE.H
index 6b3c69011c..fe57cb91f4 100644
--- a/src/ODE/ODE/ODE.H
+++ b/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;
 };
 
diff --git a/src/ODE/ODESolvers/KRR4/KRR4.C b/src/ODE/ODESolvers/KRR4/KRR4.C
index fa00919a42..ec3db04574 100644
--- a/src/ODE/ODESolvers/KRR4/KRR4.C
+++ b/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++)
         {
diff --git a/src/ODE/ODESolvers/KRR4/KRR4.H b/src/ODE/ODESolvers/KRR4/KRR4.H
index cf86ee3f15..f5393159a6 100644
--- a/src/ODE/ODESolvers/KRR4/KRR4.H
+++ b/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,
diff --git a/src/ODE/ODESolvers/SIBS/SIBS.H b/src/ODE/ODESolvers/SIBS/SIBS.H
index 75383bc38a..164b92d299 100644
--- a/src/ODE/ODESolvers/SIBS/SIBS.H
+++ b/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;
 
 
diff --git a/src/ODE/ODESolvers/SIBS/SIMPR.C b/src/ODE/ODESolvers/SIBS/SIMPR.C
index f540f26d11..671f1acec9 100644
--- a/src/ODE/ODESolvers/SIBS/SIMPR.C
+++ b/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++)
     {
diff --git a/src/ODE/ODESolvers/SIBS/polyExtrapolate.C b/src/ODE/ODESolvers/SIBS/polyExtrapolate.C
index 3f3c45d339..c22720a6b1 100644
--- a/src/ODE/ODESolvers/SIBS/polyExtrapolate.C
+++ b/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();
diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files
index 57dc7f8dc7..cf4e331e25 100644
--- a/src/OpenFOAM/Make/files
+++ b/src/OpenFOAM/Make/files
@@ -177,9 +177,9 @@ dimensionedTypes/dimensionedTensor/dimensionedTensor.C
 
 matrices/solution/solution.C
 
-scalarMatrix = matrices/scalarMatrix
-$(scalarMatrix)/scalarMatrix.C
-$(scalarMatrix)/SVD/SVD.C
+scalarMatrices = matrices/scalarMatrices
+$(scalarMatrices)/scalarMatrices.C
+$(scalarMatrices)/SVD/SVD.C
 
 LUscalarMatrix = matrices/LUscalarMatrix
 $(LUscalarMatrix)/LUscalarMatrix.C
diff --git a/src/OpenFOAM/db/IOstreams/Pstreams/PstreamCommsStruct.C b/src/OpenFOAM/db/IOstreams/Pstreams/PstreamCommsStruct.C
index 6b078aa25d..7a1a680643 100644
--- a/src/OpenFOAM/db/IOstreams/Pstreams/PstreamCommsStruct.C
+++ b/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"
diff --git a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C
index c023e59c5a..a31058eace 100644
--- a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C
+++ b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C
@@ -29,7 +29,8 @@ License
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 template<class Type>
-Foam::DiagonalMatrix<Type>::DiagonalMatrix(const Matrix<Type>& a)
+template<class Form>
+Foam::DiagonalMatrix<Type>::DiagonalMatrix(const Matrix<Form, Type>& a)
 :
     List<Type>(min(a.n(), a.m()))
 {
diff --git a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H
index 36138fc25d..ae8ceffc1c 100644
--- a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H
+++ b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H
@@ -46,7 +46,7 @@ namespace Foam
 
 // * * * * * *  * * * * * * Class Forward declaration  * * * * * * * * * * * //
 
-template<class Type> class Matrix;
+template<class Form, class Type> class Matrix;
 
 /*---------------------------------------------------------------------------*\
                            Class DiagonalMatrix Declaration
@@ -62,7 +62,8 @@ public:
     // Constructors
 
         //- Construct from diagonal component of a Matrix
-        DiagonalMatrix<Type>(const Matrix<Type>&);
+        template<class Form>
+        DiagonalMatrix<Type>(const Matrix<Form, Type>&);
 
         //- Construct empty from size
         DiagonalMatrix<Type>(const label size);
diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C
index 5de8693f2f..f425accf72 100644
--- a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C
+++ b/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, 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, 0.0);
         transfer(m);
         convert(ldum, interfaceCoeffs, interfaces);
     }
diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H
index d5e79faf19..2255cee236 100644
--- a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H
+++ b/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
diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.C b/src/OpenFOAM/matrices/Matrix/Matrix.C
index 4fbef71c2b..577cbcabcb 100644
--- a/src/OpenFOAM/matrices/Matrix/Matrix.C
+++ b/src/OpenFOAM/matrices/Matrix/Matrix.C
@@ -28,8 +28,8 @@ License
 
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
-template<class Type>
-void Foam::Matrix<Type>::allocate()
+template<class Form, class Type>
+void Foam::Matrix<Form, Type>::allocate()
 {
     if (n_ && m_)
     {
@@ -46,8 +46,8 @@ void Foam::Matrix<Type>::allocate()
 
 // * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * * //
 
-template<class Type>
-Foam::Matrix<Type>::~Matrix()
+template<class Form, class Type>
+Foam::Matrix<Form, Type>::~Matrix()
 {
     if (v_)
     {
@@ -59,16 +59,16 @@ Foam::Matrix<Type>::~Matrix()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class Type>
-const Foam::Matrix<Type>& Foam::Matrix<Type>::null()
+template<class Form, class Type>
+const Foam::Matrix<Form, Type>& Foam::Matrix<Form, Type>::null()
 {
-    Matrix<Type>* nullPtr = reinterpret_cast<Matrix<Type>*>(NULL);
+    Matrix<Form, Type>* nullPtr = reinterpret_cast<Matrix<Form, Type>*>(NULL);
     return *nullPtr;
 }
 
 
-template<class Type>
-Foam::Matrix<Type>::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<Type>::Matrix(const label n, const label m)
 {
     if (n_ < 0 || m_ < 0)
     {
-        FatalErrorIn("Matrix<Type>::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<Type>::Matrix(const label n, const label m)
 }
 
 
-template<class Type>
-Foam::Matrix<Type>::Matrix(const label n, const label m, const Type& 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<Type>::Matrix(const label n, const label m, const Type& a)
     {
         FatalErrorIn
         (
-            "Matrix<Type>::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);
     }
@@ -117,8 +117,8 @@ Foam::Matrix<Type>::Matrix(const label n, const label m, const Type& a)
 }
 
 
-template<class Type>
-Foam::Matrix<Type>::Matrix(const Matrix<Type>& a)
+template<class Form, class Type>
+Foam::Matrix<Form, Type>::Matrix(const Matrix<Form, Type>& a)
 :
     n_(a.n_),
     m_(a.m_),
@@ -139,8 +139,8 @@ Foam::Matrix<Type>::Matrix(const Matrix<Type>& a)
 }
 
 
-template<class Type>
-void Foam::Matrix<Type>::clear()
+template<class Form, class Type>
+void Foam::Matrix<Form, Type>::clear()
 {
     if (v_)
     {
@@ -153,8 +153,8 @@ void Foam::Matrix<Type>::clear()
 }
 
 
-template<class Type>
-void Foam::Matrix<Type>::transfer(Matrix<Type>& a)
+template<class Form, class Type>
+void Foam::Matrix<Form, Type>::transfer(Matrix<Form, Type>& a)
 {
     clear();
 
@@ -169,13 +169,11 @@ void Foam::Matrix<Type>::transfer(Matrix<Type>& a)
 }
 
 
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-template<class Type>
-Foam::Matrix<Type> Foam::Matrix<Type>::T() const
+template<class Form, class Type>
+Form Foam::Matrix<Form, Type>::T() const
 {
-    const Matrix<Type>& A = *this;
-    Matrix<Type> At(m(), n());
+    const Matrix<Form, Type>& A = *this;
+    Form At(m(), n());
 
     for (register label i=0; i<n(); i++)
     {
@@ -189,8 +187,10 @@ Foam::Matrix<Type> Foam::Matrix<Type>::T() const
 }
 
 
-template<class Type>
-void Foam::Matrix<Type>::operator=(const Type& t)
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+template<class Form, class Type>
+void Foam::Matrix<Form, Type>::operator=(const Type& t)
 {
     if (v_)
     {
@@ -206,12 +206,12 @@ void Foam::Matrix<Type>::operator=(const Type& t)
 
 
 // Assignment operator. Takes linear time.
-template<class Type>
-void Foam::Matrix<Type>::operator=(const Matrix<Type>& a)
+template<class Form, class Type>
+void Foam::Matrix<Form, Type>::operator=(const Matrix<Form, Type>& a)
 {
     if (this == &a)
     {
-        FatalErrorIn("Matrix<Type>::operator=(const Matrix<Type>&)")
+        FatalErrorIn("Matrix<Form, Type>::operator=(const Matrix<Form, Type>&)")
             << "attempted assignment to self"
             << abort(FatalError);
     }
@@ -240,8 +240,8 @@ void Foam::Matrix<Type>::operator=(const Matrix<Type>& a)
 
 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 
-template<class Type>
-const Type& Foam::max(const Matrix<Type>& a)
+template<class Form, class Type>
+const Type& Foam::max(const Matrix<Form, Type>& a)
 {
     label nm = a.n_*a.m_;
 
@@ -262,7 +262,7 @@ const Type& Foam::max(const Matrix<Type>& a)
     }
     else
     {
-        FatalErrorIn("max(const Matrix<Type>&)")
+        FatalErrorIn("max(const Matrix<Form, Type>&)")
             << "matrix is empty"
             << abort(FatalError);
 
@@ -272,8 +272,8 @@ const Type& Foam::max(const Matrix<Type>& a)
 }
 
 
-template<class Type>
-const Type& Foam::min(const Matrix<Type>& a)
+template<class Form, class Type>
+const Type& Foam::min(const Matrix<Form, Type>& a)
 {
     label nm = a.n_*a.m_;
 
@@ -294,7 +294,7 @@ const Type& Foam::min(const Matrix<Type>& a)
     }
     else
     {
-        FatalErrorIn("min(const Matrix<Type>&)")
+        FatalErrorIn("min(const Matrix<Form, Type>&)")
             << "matrix is empty"
             << abort(FatalError);
 
@@ -306,10 +306,10 @@ const Type& Foam::min(const Matrix<Type>& a)
 
 // * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
 
-template<class Type>
-Foam::Matrix<Type> Foam::operator-(const Matrix<Type>& a)
+template<class Form, class Type>
+Form Foam::operator-(const Matrix<Form, Type>& a)
 {
-    Matrix<Type> na(a.n_, a.m_);
+    Form na(a.n_, a.m_);
 
     if (a.n_ && a.m_)
     {
@@ -327,14 +327,14 @@ Foam::Matrix<Type> Foam::operator-(const Matrix<Type>& a)
 }
 
 
-template<class Type>
-Foam::Matrix<Type> Foam::operator+(const Matrix<Type>& a, const Matrix<Type>& 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<Type>::operator+(const Matrix<Type>&, const Matrix<Type>&)"
+            "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);
@@ -344,13 +344,13 @@ Foam::Matrix<Type> Foam::operator+(const Matrix<Type>& a, const Matrix<Type>& b)
     {
         FatalErrorIn
         (
-            "Matrix<Type>::operator+(const Matrix<Type>&, const Matrix<Type>&)"
+            "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<Type> ab(a.n_, a.m_);
+    Form ab(a.n_, a.m_);
 
     Type* abv = ab.v_[0];
     const Type* av = a.v_[0];
@@ -366,14 +366,14 @@ Foam::Matrix<Type> Foam::operator+(const Matrix<Type>& a, const Matrix<Type>& b)
 }
 
 
-template<class Type>
-Foam::Matrix<Type> Foam::operator-(const Matrix<Type>& a, const Matrix<Type>& 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<Type>::operator-(const Matrix<Type>&, const Matrix<Type>&)"
+            "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);
@@ -383,13 +383,13 @@ Foam::Matrix<Type> Foam::operator-(const Matrix<Type>& a, const Matrix<Type>& b)
     {
         FatalErrorIn
         (
-            "Matrix<Type>::operator-(const Matrix<Type>&, const Matrix<Type>&)"
+            "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<Type> ab(a.n_, a.m_);
+    Form ab(a.n_, a.m_);
 
     Type* abv = ab.v_[0];
     const Type* av = a.v_[0];
@@ -405,17 +405,17 @@ Foam::Matrix<Type> Foam::operator-(const Matrix<Type>& a, const Matrix<Type>& b)
 }
 
 
-template<class Type>
-Foam::Matrix<Type> Foam::operator*(const scalar s, const Matrix<Type>& a)
+template<class Form, class Type>
+Form Foam::operator*(const scalar s, const Matrix<Form, Type>& a)
 {
-    Matrix<Type> sa(a.n_, a.m_);
+    Form sa(a.n_, a.m_);
 
     if (a.n_ && a.m_)
     {
         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];
diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.H b/src/OpenFOAM/matrices/Matrix/Matrix.H
index 6b4da39d4f..7c63e6d530 100644
--- a/src/OpenFOAM/matrices/Matrix/Matrix.H
+++ b/src/OpenFOAM/matrices/Matrix/Matrix.H
@@ -51,37 +51,26 @@ namespace Foam
 
 // Forward declaration of friend functions and operators
 
-template<class Type> class Matrix;
+template<class Form, class Type> class Matrix;
 
-template<class Type> const Type& max(const Matrix<Type>&);
-template<class Type> const Type& min(const Matrix<Type>&);
-
-template<class Type> Matrix<Type> operator-(const Matrix<Type>&);
-template<class Type> Matrix<Type> operator+
+template<class Form, class Type> Istream& operator>>
 (
-    const Matrix<Type>&,
-    const Matrix<Type>&
-);
-template<class Type> Matrix<Type> operator-
-(
-    const Matrix<Type>&,
-    const Matrix<Type>&
-);
-template<class Type> Matrix<Type> operator*
-(
-    const scalar,
-    const Matrix<Type>&
+    Istream&,
+    Matrix<Form, Type>&
 );
 
-template<class Type> Istream& operator>>(Istream&, Matrix<Type>&);
-template<class Type> Ostream& operator<<(Ostream&, const Matrix<Type>&);
+template<class Form, class Type> Ostream& operator<<
+(
+    Ostream&,
+    const Matrix<Form, Type>&
+);
 
 
 /*---------------------------------------------------------------------------*\
                            Class Matrix Declaration
 \*---------------------------------------------------------------------------*/
 
-template<class Type>
+template<class Form, class Type>
 class Matrix
 {
     // Private data
@@ -112,13 +101,13 @@ public:
         Matrix(const label n, const label m, const Type&);
 
         //- Copy constructor.
-        Matrix(const Matrix<Type>&);
+        Matrix(const Matrix<Form, Type>&);
 
         //- Construct from Istream.
         Matrix(Istream&);
 
         //- Clone
-        inline autoPtr<Matrix<Type> > clone() const;
+        inline autoPtr<Matrix<Form, Type> > clone() const;
 
 
     // Destructor
@@ -129,7 +118,7 @@ public:
     // Member functions
 
         //- Return a null Matrix
-        static const Matrix<Type>& null();
+        static const Matrix<Form, Type>& null();
 
 
         // Access
@@ -160,7 +149,11 @@ public:
 
             //- Transfer the contents of the argument Matrix into this Matrix
             //  and annull the argument Matrix.
-            void transfer(Matrix<Type>&);
+            void transfer(Matrix<Form, Type>&);
+
+
+        //- Return the transpose of the matrix
+        Form T() const;
 
 
     // Member operators
@@ -171,52 +164,49 @@ public:
         //- Return subscript-checked element of constant Matrix.
         inline const Type* operator[](const label) const;
 
-        //- Return the transpose of the matrix
-        Matrix<Type> T() const;
-
         //- Assignment operator. Takes linear time.
-        void operator=(const Matrix<Type>&);
+        void operator=(const Matrix<Form, Type>&);
 
         //- Assignment of all entries to the given value
         void operator=(const Type&);
 
 
-    // Friend functions
-
-        friend const Type& max<Type>(const Matrix<Type>&);
-        friend const Type& min<Type>(const Matrix<Type>&);
-
-
-    // Friend operators
-
-        friend Matrix<Type> operator-<Type>(const Matrix<Type>&);
-        friend Matrix<Type> operator+<Type>
-        (
-            const Matrix<Type>&,
-            const Matrix<Type>&
-        );
-        friend Matrix<Type> operator-<Type>
-        (
-            const Matrix<Type>&,
-            const Matrix<Type>&
-        );
-        friend Matrix<Type> operator*<Type>
-        (
-            const scalar,
-            const Matrix<Type>&
-        );
-
-
     // IOstream operators
 
         //- Read Matrix from Istream, discarding contents of existing Matrix.
-        friend Istream& operator>> <Type>(Istream&, Matrix<Type>&);
+        friend Istream& operator>> <Form, Type>(Istream&, Matrix<Form, Type>&);
 
         // Write Matrix to Ostream.
-        friend Ostream& operator<< <Type>(Ostream&, const Matrix<Type>&);
+        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
diff --git a/src/OpenFOAM/matrices/Matrix/MatrixI.H b/src/OpenFOAM/matrices/Matrix/MatrixI.H
index 8323174d4f..8d758cf5b1 100644
--- a/src/OpenFOAM/matrices/Matrix/MatrixI.H
+++ b/src/OpenFOAM/matrices/Matrix/MatrixI.H
@@ -26,8 +26,8 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-template<class Type>
-inline Foam::Matrix<Type>::Matrix()
+template<class Form, class Type>
+inline Foam::Matrix<Form, Type>::Matrix()
 :
     n_(0),
     m_(0),
@@ -35,67 +35,67 @@ inline Foam::Matrix<Type>::Matrix()
 {}
 
 
-template<class Type>
-inline Foam::autoPtr<Foam::Matrix<Type> > Foam::Matrix<Type>::clone() const
+template<class Form, class Type>
+inline Foam::autoPtr<Foam::Matrix<Form, Type> > Foam::Matrix<Form, Type>::clone() const
 {
-    return autoPtr<Matrix<Type> >(new Matrix<Type>(*this));
+    return autoPtr<Matrix<Form, Type> >(new Matrix<Form, Type>(*this));
 }
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 //- Return the number of rows
-template<class Type>
-inline Foam::label Foam::Matrix<Type>::n() const
+template<class Form, class Type>
+inline Foam::label Foam::Matrix<Form, Type>::n() const
 {
     return n_;
 }
 
 
-template<class Type>
-inline Foam::label Foam::Matrix<Type>::m() const
+template<class Form, class Type>
+inline Foam::label Foam::Matrix<Form, Type>::m() const
 {
     return m_;
 }
 
 
-template<class Type>
-inline Foam::label Foam::Matrix<Type>::size() const
+template<class Form, class Type>
+inline Foam::label Foam::Matrix<Form, Type>::size() const
 {
     return n_*m_;
 }
 
 
-template<class Type>
-inline void Foam::Matrix<Type>::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<Type>::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<Type>::checki(const label)")
+        FatalErrorIn("Matrix<Form, Type>::checki(const label)")
             << "index " << i << " out of range 0 ... " << n_-1
             << abort(FatalError);
     }
 }
 
 
-template<class Type>
-inline void Foam::Matrix<Type>::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<Type>::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<Type>::checkj(const label)")
+        FatalErrorIn("Matrix<Form, Type>::checkj(const label)")
             << "index " << j << " out of range 0 ... " << m_-1
             << abort(FatalError);
     }
@@ -104,8 +104,8 @@ inline void Foam::Matrix<Type>::checkj(const label j) const
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
-template<class Type>
-inline Type* Foam::Matrix<Type>::operator[](const label i)
+template<class Form, class Type>
+inline Type* Foam::Matrix<Form, Type>::operator[](const label i)
 {
 #   ifdef FULLDEBUG
     checki(i);
@@ -114,8 +114,8 @@ inline Type* Foam::Matrix<Type>::operator[](const label i)
 }
 
 
-template<class Type>
-inline const Type* Foam::Matrix<Type>::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);
diff --git a/src/OpenFOAM/matrices/Matrix/MatrixIO.C b/src/OpenFOAM/matrices/Matrix/MatrixIO.C
index 9b5a922f0d..15f0d5345e 100644
--- a/src/OpenFOAM/matrices/Matrix/MatrixIO.C
+++ b/src/OpenFOAM/matrices/Matrix/MatrixIO.C
@@ -32,8 +32,8 @@ License
 
 // * * * * * * * * * * * * * * * Ostream Operator *  * * * * * * * * * * * * //
 
-template<class Type>
-Foam::Matrix<Type>::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<Type>::Matrix(Istream& is)
 }
 
 
-template<class Type>
-Foam::Istream& Foam::operator>>(Istream& is, Matrix<Type>& 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<Type>&)");
+    is.fatalCheck("operator>>(Istream&, Matrix<Form, Type>&)");
 
     token firstToken(is);
 
-    is.fatalCheck("operator>>(Istream&, Matrix<Type>&) : reading first token");
+    is.fatalCheck("operator>>(Istream&, Matrix<Form, Type>&) : reading first token");
 
     if (firstToken.isLabel())
     {
@@ -88,7 +88,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<Type>& M)
 
                             is.fatalCheck
                             (
-                                "operator>>(Istream&, Matrix<Type>&) : "
+                                "operator>>(Istream&, Matrix<Form, Type>&) : "
                                 "reading entry"
                             );
                         }
@@ -103,7 +103,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<Type>& M)
 
                     is.fatalCheck
                     (
-                        "operator>>(Istream&, Matrix<Type>&) : "
+                        "operator>>(Istream&, Matrix<Form, Type>&) : "
                         "reading the single entry"
                     );
 
@@ -128,7 +128,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<Type>& M)
 
                 is.fatalCheck
                 (
-                    "operator>>(Istream&, Matrix<Type>&) : "
+                    "operator>>(Istream&, Matrix<Form, Type>&) : "
                     "reading the binary block"
                 );
             }
@@ -136,7 +136,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<Type>& M)
     }
     else
     {
-        FatalIOErrorIn("operator>>(Istream&, Matrix<Type>&)", is)
+        FatalIOErrorIn("operator>>(Istream&, Matrix<Form, Type>&)", is)
             << "incorrect first token, expected <int>, found "
             << firstToken.info()
             << exit(FatalIOError);
@@ -146,8 +146,8 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<Type>& M)
 }
 
 
-template<class Type>
-Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<Type>& M)
+template<class Form, class Type>
+Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<Form, Type>& M)
 {
     label nm = M.n_*M.m_;
 
diff --git a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H
new file mode 100644
index 0000000000..e5e6b16777
--- /dev/null
+++ b/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
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H
new file mode 100644
index 0000000000..1cacd2f11a
--- /dev/null
+++ b/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));
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H
new file mode 100644
index 0000000000..6239781d85
--- /dev/null
+++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H
@@ -0,0 +1,96 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::SquareMatrix
+
+Description
+    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
+    SquareMatrixI.H
+    SquareMatrix.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef SquareMatrix_H
+#define SquareMatrix_H
+
+#include "Matrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class Matrix Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class SquareMatrix
+:
+    public Matrix<SquareMatrix<Type>, Type>
+{
+
+public:
+
+    // Constructors
+
+        //- Null constructor.
+        inline SquareMatrix();
+
+        //- Construct given number of rows/columns.
+        inline SquareMatrix(const label n);
+
+        //- Construct given number of rows and columns,
+        //  It checks that m == n.
+        inline SquareMatrix(const label m, const label n);
+
+        //- Construct with given number of rows/columns
+        //  and value for all elements.
+        inline SquareMatrix(const label n, const Type&);
+
+        //- Construct from Istream.
+        inline SquareMatrix(Istream&);
+
+        //- Clone
+        inline autoPtr<SquareMatrix<Type> > clone() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#   include "SquareMatrixI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H
new file mode 100644
index 0000000000..a7371a6a2e
--- /dev/null
+++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H
@@ -0,0 +1,75 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 n, const Type& t)
+:
+    Matrix<SquareMatrix<Type>, Type>(n, t)
+{}
+
+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));
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.C b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C
similarity index 96%
rename from src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.C
rename to src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C
index a1355c393d..13ef29ed9b 100644
--- a/src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.C
+++ b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C
@@ -26,12 +26,12 @@ License
 
 #include "SVD.H"
 #include "scalarList.H"
-#include "scalarMatrix.H"
+#include "scalarMatrices.H"
 #include "ListOps.H"
 
 // * * * * * * * * * * * * * * * * Constructor  * * * * * * * * * * * * * * //
 
-Foam::SVD::SVD(const Matrix<scalar>& A, const scalar minCondition)
+Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
 :
     U_(A),
     V_(A.m(), A.m()),
@@ -294,8 +294,11 @@ Foam::SVD::SVD(const Matrix<scalar>& A, const scalar minCondition)
             }
             if (its == 34)
             {
-                WarningIn("SVD::SVD(Matrix<scalar>& A)")
-                    << "no convergence in 35 SVD iterations"
+                WarningIn
+                (
+                    "SVD::SVD"
+                    "(scalarRectangularMatrix& A, const scalar minCondition)"
+                )   << "no convergence in 35 SVD iterations"
                     << endl;
             }
 
@@ -375,7 +378,7 @@ Foam::SVD::SVD(const Matrix<scalar>& A, const scalar minCondition)
     multiply(VSinvUt_, V_, inv(S_), U_.T());
 
     // test SVD
-    /*Matrix<scalar> SVDA(A.n(), A.m());
+    /*scalarRectangularMatrix SVDA(A.n(), A.m());
     multiply(SVDA, U_, S_, transpose(V_));
     scalar maxDiff = 0;
     scalar diff = 0;
diff --git a/src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.H b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.H
similarity index 86%
rename from src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.H
rename to src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.H
index b597027732..48e581198a 100644
--- a/src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.H
+++ b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.H
@@ -23,12 +23,13 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    SVD
+    Foam::SVD
 
 Description
     Singular value decomposition of a rectangular matrix.
 
 SourceFiles
+    SVDI.H
     SVD.C
 
 \*---------------------------------------------------------------------------*/
@@ -36,8 +37,7 @@ SourceFiles
 #ifndef SVD_H
 #define SVD_H
 
-#include "DiagonalMatrix.H"
-#include "Matrix.H"
+#include "scalarMatrices.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -56,16 +56,16 @@ class SVD
     // Private data
 
         //- Rectangular matrix with the same dimensions as the input
-        Matrix<scalar> U_;
+        scalarRectangularMatrix U_;
 
         //- square matrix V
-        Matrix<scalar> V_;
+        scalarRectangularMatrix V_;
 
         //- The singular values
         DiagonalMatrix<scalar> S_;
 
         //- The matrix product V S^(-1) U^T
-        Matrix<scalar> VSinvUt_;
+        scalarRectangularMatrix VSinvUt_;
 
         //- The number of zero singular values
         label nZeros_;
@@ -88,22 +88,22 @@ public:
     // Constructors
 
         //- Construct from a rectangular Matrix
-        SVD(const Matrix<scalar>& A, const scalar minCondition = 0);
+        SVD(const scalarRectangularMatrix& A, const scalar minCondition = 0);
 
 
     // Access functions
 
         //- Return U
-        inline const Matrix<scalar>& U() const;
+        inline const scalarRectangularMatrix& U() const;
 
         //- Return the square matrix V
-        inline const Matrix<scalar>& V() const;
+        inline const scalarRectangularMatrix& V() const;
 
         //- Return the singular values
-        inline const DiagonalMatrix<scalar>& S() const;
+        inline const scalarDiagonalMatrix& S() const;
 
         //- Return VSinvUt (the pseudo inverse)
-        inline const Matrix<scalar>& VSinvUt() const;
+        inline const scalarRectangularMatrix& VSinvUt() const;
 
         //- Return the number of zero singular values
         inline label nZeros() const;
diff --git a/src/OpenFOAM/matrices/scalarMatrix/SVD/SVDI.H b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVDI.H
similarity index 88%
rename from src/OpenFOAM/matrices/scalarMatrix/SVD/SVDI.H
rename to src/OpenFOAM/matrices/scalarMatrices/SVD/SVDI.H
index 337647edfb..298aac1aeb 100644
--- a/src/OpenFOAM/matrices/scalarMatrix/SVD/SVDI.H
+++ b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVDI.H
@@ -35,22 +35,22 @@ inline const T Foam::SVD::sign(const T& a, const T& b)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline const Foam::Matrix<Foam::scalar>& Foam::SVD::U() const
+inline const Foam::scalarRectangularMatrix& Foam::SVD::U() const
 {
     return U_;
 }
 
-inline const Foam::Matrix<Foam::scalar>& Foam::SVD::V() const
+inline const Foam::scalarRectangularMatrix& Foam::SVD::V() const
 {
     return V_;
 }
 
-inline const Foam::DiagonalMatrix<Foam::scalar>& Foam::SVD::S() const
+inline const Foam::scalarDiagonalMatrix& Foam::SVD::S() const
 {
     return S_;
 }
 
-inline const Foam::Matrix<Foam::scalar>& Foam::SVD::VSinvUt() const
+inline const Foam::scalarRectangularMatrix& Foam::SVD::VSinvUt() const
 {
     return VSinvUt_;
 }
diff --git a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.C b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C
similarity index 76%
rename from src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.C
rename to src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C
index b88e167b6f..8ddd01a82c 100644
--- a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.C
+++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C
@@ -24,39 +24,14 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "scalarMatrix.H"
+#include "scalarMatrices.H"
 #include "SVD.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
+void Foam::LUDecompose
 (
-    Matrix<scalar>& matrix,
+    scalarSquareMatrix& matrix,
     labelList& pivotIndices
 )
 {
@@ -81,8 +56,8 @@ void Foam::scalarMatrix::LUDecompose
         {
             FatalErrorIn
             (
-                "scalarMatrix::LUdecompose"
-                "(Matrix<scalar>& matrix, labelList& rowIndices)"
+                "LUdecompose"
+                "(scalarSquareMatrix& matrix, labelList& rowIndices)"
             )   << "Singular matrix" << exit(FatalError);
         }
 
@@ -164,9 +139,9 @@ void Foam::scalarMatrix::LUDecompose
 
 void Foam::multiply
 (
-    Matrix<scalar>& ans,         // value changed in return
-    const Matrix<scalar>& A,
-    const Matrix<scalar>& B
+    scalarRectangularMatrix& ans,         // value changed in return
+    const scalarRectangularMatrix& A,
+    const scalarRectangularMatrix& B
 )
 {
     if (A.m() != B.n())
@@ -174,15 +149,15 @@ void Foam::multiply
         FatalErrorIn
         (
             "multiply("
-            "Matrix<scalar>& answer "
-            "const Matrix<scalar>& A, "
-            "const Matrix<scalar>& B)"
+            "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 = Matrix<scalar>(A.n(), B.m(), scalar(0));
+    ans = scalarRectangularMatrix(A.n(), B.m(), scalar(0));
 
     for(register label i = 0; i < A.n(); i++)
     {
@@ -199,10 +174,10 @@ void Foam::multiply
 
 void Foam::multiply
 (
-    Matrix<scalar>& ans,         // value changed in return
-    const Matrix<scalar>& A,
-    const Matrix<scalar>& B,
-    const Matrix<scalar>& C
+    scalarRectangularMatrix& ans,         // value changed in return
+    const scalarRectangularMatrix& A,
+    const scalarRectangularMatrix& B,
+    const scalarRectangularMatrix& C
 )
 {
     if (A.m() != B.n())
@@ -210,10 +185,10 @@ void Foam::multiply
         FatalErrorIn
         (
             "multiply("
-            "const Matrix<scalar>& A, "
-            "const Matrix<scalar>& B, "
-            "const Matrix<scalar>& C, "
-            "Matrix<scalar>& answer)"
+            "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);
@@ -224,16 +199,16 @@ void Foam::multiply
         FatalErrorIn
         (
             "multiply("
-            "const Matrix<scalar>& A, "
-            "const Matrix<scalar>& B, "
-            "const Matrix<scalar>& C, "
-            "Matrix<scalar>& answer)"
+            "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 = Matrix<scalar>(A.n(), C.m(), scalar(0));
+    ans = scalarRectangularMatrix(A.n(), C.m(), scalar(0));
 
     for(register label i = 0; i < A.n(); i++)
     {
@@ -255,10 +230,10 @@ void Foam::multiply
 
 void Foam::multiply
 (
-    Matrix<scalar>& ans,         // value changed in return
-    const Matrix<scalar>& A,
+    scalarRectangularMatrix& ans,         // value changed in return
+    const scalarRectangularMatrix& A,
     const DiagonalMatrix<scalar>& B,
-    const Matrix<scalar>& C
+    const scalarRectangularMatrix& C
 )
 {
     if (A.m() != B.size())
@@ -266,10 +241,10 @@ void Foam::multiply
         FatalErrorIn
         (
             "multiply("
-            "const Matrix<scalar>& A, "
+            "const scalarRectangularMatrix& A, "
             "const DiagonalMatrix<scalar>& B, "
-            "const Matrix<scalar>& C, "
-            "Matrix<scalar>& answer)"
+            "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);
@@ -280,16 +255,16 @@ void Foam::multiply
         FatalErrorIn
         (
             "multiply("
-            "const Matrix<scalar>& A, "
+            "const scalarRectangularMatrix& A, "
             "const DiagonalMatrix<scalar>& B, "
-            "const Matrix<scalar>& C, "
-            "Matrix<scalar>& answer)"
+            "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 = Matrix<scalar>(A.n(), C.m(), scalar(0));
+    ans = scalarRectangularMatrix(A.n(), C.m(), scalar(0));
 
     for(register label i = 0; i < A.n(); i++)
     {
@@ -304,9 +279,9 @@ void Foam::multiply
 }
 
 
-Foam::Matrix<Foam::scalar> Foam::SVDinv
+Foam::RectangularMatrix<Foam::scalar> Foam::SVDinv
 (
-    const Matrix<scalar>& A,
+    const scalarRectangularMatrix& A,
     scalar minCondition
 )
 {
diff --git a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H
new file mode 100644
index 0000000000..17d8727516
--- /dev/null
+++ b/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
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrixTemplates.C b/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C
similarity index 90%
rename from src/OpenFOAM/matrices/scalarMatrix/scalarMatrixTemplates.C
rename to src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C
index fb48a9b26d..131856b3fc 100644
--- a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrixTemplates.C
+++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C
@@ -24,15 +24,15 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "scalarMatrix.H"
+#include "scalarMatrices.H"
 #include "Swap.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
-void Foam::scalarMatrix::solve
+void Foam::solve
 (
-    Matrix<scalar>& tmpMatrix,
+    scalarSquareMatrix& tmpMatrix,
     Field<Type>& sourceSol
 )
 {
@@ -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);
         }
@@ -102,18 +102,23 @@ void Foam::scalarMatrix::solve
 
 
 template<class Type>
-void Foam::scalarMatrix::solve(Field<Type>& psi, const Field<Type>& source) const
+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 Type>
-void Foam::scalarMatrix::LUBacksubstitute
+void Foam::LUBacksubstitute
 (
-    const Matrix<scalar>& luMatrix,
+    const scalarSquareMatrix& luMatrix,
     const labelList& pivotIndices,
     Field<Type>& sourceSol
 )
@@ -160,9 +165,9 @@ void Foam::scalarMatrix::LUBacksubstitute
 
 
 template<class Type>
-void Foam::scalarMatrix::LUsolve
+void Foam::LUsolve
 (
-    Matrix<scalar>& matrix,
+    scalarSquareMatrix& matrix,
     Field<Type>& sourceSol
 )
 {
diff --git a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.H b/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.H
deleted file mode 100644
index 2d22e71d4e..0000000000
--- a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.H
+++ /dev/null
@@ -1,155 +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
-
-Class
-    Foam::scalarMatrix
-
-Description
-    Foam::scalarMatrix
-
-SourceFiles
-    scalarMatrix.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef scalarMatrix_H
-#define scalarMatrix_H
-
-#include "Matrix.H"
-#include "DiagonalMatrix.H"
-#include "scalarField.H"
-#include "labelList.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                           Class scalarMatrix Declaration
-\*---------------------------------------------------------------------------*/
-
-class scalarMatrix
-:
-    public Matrix<scalar>
-{
-
-public:
-
-    // Constructors
-
-        //- Construct null
-        scalarMatrix();
-
-        //- Construct given size
-        scalarMatrix(const label);
-
-        //- Construct from Matrix<scalar>
-        scalarMatrix(const Matrix<scalar>&);
-
-        //- Construct from Istream
-        scalarMatrix(Istream&);
-
-
-    // Member Functions
-
-        //- Solve the matrix using Gaussian elimination with pivoting,
-        //  returning the solution in the source
-        template<class Type>
-        static void solve(Matrix<scalar>& 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 Field<Type>& 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 Type>
-        static void LUBacksubstitute
-        (
-            const Matrix<scalar>& 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>
-        static void LUsolve(Matrix<scalar>& matrix, Field<Type>& source);
-};
-
-
-// Global functions
-
-void multiply
-(
-    Matrix<scalar>& answer,         // value changed in return
-    const Matrix<scalar>& A,
-    const Matrix<scalar>& B
-);
-
-void multiply
-(
-    Matrix<scalar>& answer,         // value changed in return
-    const Matrix<scalar>& A,
-    const Matrix<scalar>& B,
-    const Matrix<scalar>& C
-);
-
-void multiply
-(
-    Matrix<scalar>& answer,         // value changed in return
-    const Matrix<scalar>& A,
-    const DiagonalMatrix<scalar>& B,
-    const Matrix<scalar>& C
-);
-
-//- Return the inverse of matrix A using SVD
-Matrix<scalar> SVDinv(const Matrix<scalar>& A, scalar minCondition = 0);
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-#   include "scalarMatrixTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C
index 1265b87ee5..f0f80750a5 100644
--- a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C
+++ b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C
@@ -31,7 +31,7 @@ License
 template<class Type>
 Foam::simpleMatrix<Type>::simpleMatrix(const label mSize)
 :
-    scalarMatrix(mSize),
+    scalarSquareMatrix(mSize),
     source_(mSize, pTraits<Type>::zero)
 {}
 
@@ -39,11 +39,11 @@ Foam::simpleMatrix<Type>::simpleMatrix(const label mSize)
 template<class Type>
 Foam::simpleMatrix<Type>::simpleMatrix
 (
-    const scalarMatrix& matrix,
+    const scalarSquareMatrix& matrix,
     const Field<Type>& source
 )
 :
-    scalarMatrix(matrix),
+    scalarSquareMatrix(matrix),
     source_(source)
 {}
 
@@ -51,7 +51,7 @@ Foam::simpleMatrix<Type>::simpleMatrix
 template<class Type>
 Foam::simpleMatrix<Type>::simpleMatrix(Istream& is)
 :
-    scalarMatrix(is),
+    scalarSquareMatrix(is),
     source_(is)
 {}
 
@@ -61,10 +61,10 @@ Foam::simpleMatrix<Type>::simpleMatrix(Istream& is)
 template<class Type>
 Foam::Field<Type> Foam::simpleMatrix<Type>::solve() const
 {
-    scalarMatrix tmpMatrix = *this;
+    scalarSquareMatrix tmpMatrix = *this;
     Field<Type> sourceSol = source_;
 
-    scalarMatrix::solve(tmpMatrix, sourceSol);
+    Foam::solve(tmpMatrix, sourceSol);
 
     return sourceSol;
 }
@@ -73,10 +73,10 @@ Foam::Field<Type> Foam::simpleMatrix<Type>::solve() const
 template<class Type>
 Foam::Field<Type> Foam::simpleMatrix<Type>::LUsolve() const
 {
-    scalarMatrix luMatrix = *this;
+    scalarSquareMatrix luMatrix = *this;
     Field<Type> sourceSol = source_;
 
-    scalarMatrix::LUsolve(luMatrix, sourceSol);
+    Foam::LUsolve(luMatrix, sourceSol);
 
     return sourceSol;
 }
@@ -108,7 +108,7 @@ void Foam::simpleMatrix<Type>::operator=(const simpleMatrix<Type>& m)
             << abort(FatalError);
     }
 
-    scalarMatrix::operator=(m);
+    scalarSquareMatrix::operator=(m);
     source_ = m.source_;
 }
 
@@ -124,8 +124,8 @@ Foam::simpleMatrix<Type> Foam::operator+
 {
     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_
     );
 }
@@ -140,8 +140,8 @@ Foam::simpleMatrix<Type> Foam::operator-
 {
     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_
     );
 }
@@ -159,7 +159,7 @@ Foam::simpleMatrix<Type> Foam::operator*(const scalar s, const simpleMatrix<Type
 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;
 }
 
diff --git a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H
index d7b1a027be..a972c787c8 100644
--- a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H
+++ b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H
@@ -36,7 +36,7 @@ SourceFiles
 #ifndef simpleMatrix_H
 #define simpleMatrix_H
 
-#include "scalarMatrix.H"
+#include "scalarMatrices.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -48,27 +48,6 @@ namespace Foam
 template<class Type>
 class simpleMatrix;
 
-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>&
-);
-
 template<class Type>
 Ostream& operator<<
 (
@@ -84,7 +63,7 @@ Ostream& operator<<
 template<class Type>
 class simpleMatrix
 :
-    public scalarMatrix
+    public scalarSquareMatrix
 {
     // Private data
 
@@ -99,7 +78,7 @@ public:
         simpleMatrix(const label);
 
         //- Construct from components
-        simpleMatrix(const scalarMatrix&, const Field<Type>&);
+        simpleMatrix(const scalarSquareMatrix&, const Field<Type>&);
 
         //- Construct from Istream
         simpleMatrix(Istream&);
@@ -137,27 +116,6 @@ public:
         void operator=(const simpleMatrix<Type>&);
 
 
-    // Friend Operators
-
-        friend simpleMatrix<Type> operator+ <Type>
-        (
-            const simpleMatrix<Type>&,
-            const simpleMatrix<Type>&
-        );
-
-        friend simpleMatrix<Type> operator- <Type>
-        (
-            const simpleMatrix<Type>&,
-            const simpleMatrix<Type>&
-        );
-
-        friend simpleMatrix<Type> operator* <Type>
-        (
-            const scalar,
-            const simpleMatrix<Type>&
-        );
-
-
     // Ostream Operator
 
         friend Ostream& operator<< <Type>
@@ -168,6 +126,30 @@ public:
 };
 
 
+// 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
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C
index fd0e55fcdd..0caebcd421 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C
@@ -49,11 +49,11 @@ Foam::quadraticFitData::quadraticFitData
 :
     MeshObject<fvMesh, quadraticFitData>(mesh),
     centralWeight_(cWeight),
-#ifdef SPHERICAL_GEOMETRY
+#   ifdef SPHERICAL_GEOMETRY
     dim_(2),
-#else
+#   else
     dim_(mesh.nGeometricD()),
-#endif
+#   endif
     minSize_
     (
         dim_ == 1 ? 3 :
@@ -204,7 +204,7 @@ Foam::label Foam::quadraticFitData::calcFit
     scalar scale = 0;
 
     // calculate the matrix of the polynomial components
-    Matrix<scalar> B(C.size(), minSize_, scalar(0));
+    scalarRectangularMatrix B(C.size(), minSize_, scalar(0));
 
     for(label ip = 0; ip < C.size(); ip++)
     {
@@ -212,11 +212,11 @@ Foam::label Foam::quadraticFitData::calcFit
 
         scalar px = (p - p0)&idir;
         scalar py = (p - p0)&jdir;
-#ifndef SPHERICAL_GEOMETRY
+#       ifndef SPHERICAL_GEOMETRY
         scalar pz = (p - p0)&kdir;
-#else
+#       else
         scalar pz = mag(p) - mag(p0);
-#endif
+#       endif
 
         if (ip == 0)
         {
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.C
index 3dd49e126e..f2da51c746 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.C
+++ b/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++)
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.H
index 512f606e5f..354dd60b4f 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.H
+++ b/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