Merge branch 'master' into cvm

applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C

@@ -462,6 +462,7 @@ int main(int argc, char *argv[])
 
 #   include "setRootCase.H"
 #   include "createTime.H"
+    runTime.functionObjects().off();
 #   include "createPolyMesh.H"
 
     scalar minLen(readScalar(IStringStream(args.additionalArgs()[0])()));

applications/utilities/mesh/advanced/combinePatchFaces/combinePatchFaces.C

@@ -439,6 +439,7 @@ int main(int argc, char *argv[])
 
 #   include "setRootCase.H"
 #   include "createTime.H"
+    runTime.functionObjects().off();
 #   include "createPolyMesh.H"
 
     scalar featureAngle(readScalar(IStringStream(args.additionalArgs()[0])()));

applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C

@@ -332,6 +332,7 @@ int main(int argc, char *argv[])
 
 #   include "setRootCase.H"
 #   include "createTime.H"
+    runTime.functionObjects().off();
 #   include "createPolyMesh.H"
 
     bool overwrite = args.options().found("overwrite");

applications/utilities/mesh/advanced/refineHexMesh/refineHexMesh.C

@@ -56,6 +56,7 @@ int main(int argc, char *argv[])
     argList::validArgs.append("cellSet");
 #   include "setRootCase.H"
 #   include "createTime.H"
+    runTime.functionObjects().off();
 #   include "createMesh.H"
     pointMesh pMesh(mesh);
 

applications/utilities/mesh/advanced/refineWallLayer/refineWallLayer.C

@@ -54,6 +54,7 @@ int main(int argc, char *argv[])
 
 #   include "setRootCase.H"
 #   include "createTime.H"
+    runTime.functionObjects().off();
 #   include "createPolyMesh.H"
 
     word patchName(args.additionalArgs()[0]);

applications/utilities/mesh/advanced/removeFaces/removeFaces.C

@@ -53,6 +53,7 @@ int main(int argc, char *argv[])
 
 #   include "setRootCase.H"
 #   include "createTime.H"
+    runTime.functionObjects().off();
 #   include "createMesh.H"
 
     bool overwrite = args.options().found("overwrite");

applications/utilities/mesh/advanced/splitCells/splitCells.C

@@ -532,6 +532,7 @@ int main(int argc, char *argv[])
 
 #   include "setRootCase.H"
 #   include "createTime.H"
+    runTime.functionObjects().off();
 #   include "createPolyMesh.H"
 
     scalar featureAngle(readScalar(IStringStream(args.additionalArgs()[0])()));

applications/utilities/mesh/conversion/polyDualMesh/polyDualMeshApp.C

@@ -47,6 +47,7 @@ int main(int argc, char *argv[])
 
 #   include "setRootCase.H"
 #   include "createTime.H"
+    runTime.functionObjects().off();
 #   include "createPolyMesh.H"
 
     scalar featureAngle(readScalar(IStringStream(args.additionalArgs()[0])()));

applications/utilities/mesh/conversion/writeMeshObj/writeMeshObj.C

@@ -346,6 +346,7 @@ int main(int argc, char *argv[])
 #   include "addTimeOptions.H"
 #   include "setRootCase.H"
 #   include "createTime.H"
+    runTime.functionObjects().off();
 
     bool patchFaces = args.options().found("patchFaces");
     bool doCell = args.options().found("cell");

applications/utilities/mesh/generation/extrude2DMesh/doExtrude2DMesh.C

@@ -61,6 +61,7 @@ int main(int argc, char *argv[])
     argList::validOptions.insert("overwrite", "");
 #   include "setRootCase.H"
 #   include "createTime.H"
+    runTime.functionObjects().off();
 #   include "createPolyMesh.H"
 
     scalar thickness(readScalar(IStringStream(args.additionalArgs()[0])()));

applications/utilities/mesh/manipulation/attachMesh/attachMesh.C

@@ -46,6 +46,7 @@ int main(int argc, char *argv[])
 
 #   include "setRootCase.H"
 #   include "createTime.H"
+    runTime.functionObjects().off();
 #   include "createPolyMesh.H"
 
     bool overwrite = args.options().found("overwrite");

applications/utilities/mesh/manipulation/autoPatch/autoPatch.C

@@ -75,6 +75,7 @@ int main(int argc, char *argv[])
 
 #   include "setRootCase.H"
 #   include "createTime.H"
+    runTime.functionObjects().off();
 #   include "createPolyMesh.H"
 
     Info<< "Mesh read in = "

applications/utilities/mesh/manipulation/createBaffles/createBaffles.C

@@ -58,6 +58,7 @@ int main(int argc, char *argv[])
 
 #   include "setRootCase.H"
 #   include "createTime.H"
+    runTime.functionObjects().off();
 #   include "createMesh.H"
 
     const polyBoundaryMesh& patches = mesh.boundaryMesh();

applications/utilities/mesh/manipulation/createPatch/createPatch.C

@@ -308,6 +308,7 @@ int main(int argc, char *argv[])
 
 #   include "setRootCase.H"
 #   include "createTime.H"
+    runTime.functionObjects().off();
 
     const bool overwrite = args.options().found("overwrite");
 

applications/utilities/mesh/manipulation/mergeOrSplitBaffles/mergeOrSplitBaffles.C

@@ -160,6 +160,7 @@ int main(int argc, char *argv[])
     argList::validOptions.insert("overwrite", "");
 #   include "setRootCase.H"
 #   include "createTime.H"
+    runTime.functionObjects().off();
 #   include "createMesh.H"
 
     bool split = args.options().found("split");

applications/utilities/mesh/manipulation/refineMesh/refineMesh.C

@@ -298,6 +298,7 @@ int main(int argc, char *argv[])
 
 #   include "setRootCase.H"
 #   include "createTime.H"
+    runTime.functionObjects().off();
 #   include "createPolyMesh.H"
 
     printEdgeStats(mesh);

applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C

@@ -374,6 +374,7 @@ int main(int argc, char *argv[])
 
 #   include "setRootCase.H"
 #   include "createTime.H"
+    runTime.functionObjects().off();
 
     // Get times list
     instantList Times = runTime.times();

applications/utilities/mesh/manipulation/splitMesh/splitMesh.C

@@ -122,6 +122,7 @@ int main(int argc, char *argv[])
 
 #   include "setRootCase.H"
 #   include "createTime.H"
+    runTime.functionObjects().off();
 #   include "createPolyMesh.H"
 
     word setName(args.additionalArgs()[0]);

applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C

@@ -1122,6 +1122,7 @@ int main(int argc, char *argv[])
 
 #   include "setRootCase.H"
 #   include "createTime.H"
+    runTime.functionObjects().off();
 #   include "createMesh.H"
 
     word blockedFacesName;

applications/utilities/mesh/manipulation/stitchMesh/stitchMesh.C

@@ -135,6 +135,7 @@ int main(int argc, char *argv[])
 
 #   include "setRootCase.H"
 #   include "createTime.H"
+    runTime.functionObjects().off();
 #   include "createMesh.H"
 
 

applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C

@@ -157,6 +157,7 @@ int main(int argc, char *argv[])
 
 #   include "setRootCase.H"
 #   include "createTime.H"
+    runTime.functionObjects().off();
 #   include "createMesh.H"
 
     word setName(args.additionalArgs()[0]);

applications/utilities/miscellaneous/foamInfoExec/foamInfoExec.C

@@ -100,13 +100,23 @@ int main(int argc, char *argv[])
 
                 if (dict.found(entryNames[0]))
                 {
-                    const entry* entPtr = &dict.lookupEntry(entryNames[0]);
+                    const entry* entPtr = &dict.lookupEntry
+                    (
+                        entryNames[0],
+                        false,
+                        true            // wildcards
+                    );
 
                     for (int i=1; i<entryNames.size(); i++)
                     {
                         if (entPtr->dict().found(entryNames[i]))
                         {
-                            entPtr = &entPtr->dict().lookupEntry(entryNames[i]);
+                            entPtr = &entPtr->dict().lookupEntry
+                            (
+                                entryNames[i],
+                                false,
+                                true    // wildcards
+                            );
                         }
                         else
                         {

applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructor.H

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

applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructorReconstructFields.C

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

applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C

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

applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamI.H

@@ -37,7 +37,7 @@ inline Foam::word Foam::vtkPV3Foam::getFirstWord(const char* str)
         {
             ++n;
         }
-        return word(str, n);
+        return word(str, n, true);
     }
     else
     {

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

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

applications/utilities/postProcessing/miscellaneous/postChannel/postChannelDict

@@ -22,6 +22,7 @@ Ny
  25 
 ); 
 Nz 30; 
- 
+
+symmetric true;
 
 // ************************************************************************* //

applications/utilities/preProcessing/changeDictionary/changeDictionary.C

@@ -164,7 +164,16 @@ int main(int argc, char *argv[])
             forAll(dictList, i)
             {
                 doneKeys[i] = dictList[i].keyword();
-                dictList.set(i, fieldDict.lookupEntry(doneKeys[i]).clone());
+                dictList.set
+                (
+                    i,
+                    fieldDict.lookupEntry
+                    (
+                        doneKeys[i],
+                        false,
+                        true
+                    ).clone()
+                );
                 fieldDict.remove(doneKeys[i]);
             }
             // Add remaining entries

applications/utilities/surface/surfaceSubset/surfaceSubset.C

@@ -100,6 +100,11 @@ int main(int argc, char *argv[])
         meshSubsetDict.lookup("addFaceNeighbours")
     );
 
+    Switch invertSelection
+    (
+        meshSubsetDict.lookup("invertSelection")
+    );
+
     // Mark the cells for the subset
 
     // Faces to subset
@@ -337,6 +342,27 @@ int main(int argc, char *argv[])
             << " faces because of addFaceNeighbours" << endl;
     }
 
+
+    if (invertSelection)
+    {
+        Info<< "Inverting selection." << endl;
+        boolList newFacesToSubset(facesToSubset.size());
+
+        forAll(facesToSubset, i)
+        {
+            if (facesToSubset[i])
+            {
+                newFacesToSubset[i] = false;
+            }
+            else
+            {
+                newFacesToSubset[i] = true;
+            }
+        }
+        facesToSubset.transfer(newFacesToSubset);
+    }
+
+
     // Create subsetted surface
     labelList pointMap;
     labelList faceMap;

bin/mpirunDebug

@@ -70,15 +70,17 @@ if [ ! -x "$exec" ]; then
     exit 1
 fi
 
+if [ ! "$PWD" ]; then
+    PWD=`pwd`
+fi
 
+echo "run $args" > $PWD/gdbCommands
+echo "where" >> $PWD/gdbCommands
+echo "Constructed gdb initialization file $PWD/gdbCommands"
 
-echo "run $args" > $HOME/gdbCommands
+$ECHO "Choose running method: 0)normal  1)gdb+xterm  2)gdb  3)log  4)log+xterm  5)xterm+valgrind: \c"
-echo "where" >> $HOME/gdbCommands
-echo "Constructed gdb initialization file $HOME/gdbCommands"
-
-$ECHO "Choose running method: 1)gdb+xterm  2)gdb  3)log  4)log+xterm  5)xterm+valgrind: \c"
 read method
-if [ "$method" -ne 1 -a "$method" -ne 2 -a "$method" -ne 3 -a "$method" -ne 4 -a "$method" -ne 5 ]; then
+if [ "$method" -ne 0 -a "$method" -ne 1 -a "$method" -ne 2 -a "$method" -ne 3 -a "$method" -ne 4 -a "$method" -ne 5 ]; then
     printUsage
     exit 1
 fi
@@ -119,15 +121,15 @@ fi
 
 echo "**sourceFoam:$sourceFoam"
 
-rm -f $HOME/mpirun.schema
+rm -f $PWD/mpirun.schema
-touch $HOME/mpirun.schema
+touch $PWD/mpirun.schema
 
 proc=0
 xpos=0
 ypos=0
 for ((proc=0; proc<$nProcs; proc++))
 do
-    procCmdFile="$HOME/processor${proc}.sh"
+    procCmdFile="$PWD/processor${proc}.sh"
     procLog="processor${proc}.log"
     geom="-geometry 120x20+$xpos+$ypos"
     node=""
@@ -141,22 +143,25 @@ do
     fi
 
     echo "#!/bin/sh" > $procCmdFile
-    if [ "$method" -eq 1 ]; then
+    if [ "$method" -eq 0 ]; then
-        echo "$sourceFoam; cd $PWD; gdb -command $HOME/gdbCommands $exec 2>&1 | tee $procLog; read dummy" >> $procCmdFile
+        echo "$sourceFoam; cd $PWD; $exec $args | tee $procLog" >> $procCmdFile
+        echo "${node}$procCmdFile" >> $PWD/mpirun.schema
+    elif [ "$method" -eq 1 ]; then
+        echo "$sourceFoam; cd $PWD; gdb -command $PWD/gdbCommands $exec 2>&1 | tee $procLog; read dummy" >> $procCmdFile
         #echo "$sourceFoam; cd $PWD; $exec $args; read dummy" >> $procCmdFile
-        echo "${node}xterm -font fixed -title 'processor'$proc $geom -e $procCmdFile" >> $HOME/mpirun.schema
+        echo "${node}xterm -font fixed -title 'processor'$proc $geom -e $procCmdFile" >> $PWD/mpirun.schema
     elif [ "$method" -eq 2 ]; then
-        echo "$sourceFoam; cd $PWD; gdb -command $HOME/gdbCommands >& $procLog" >> $procCmdFile
+        echo "$sourceFoam; cd $PWD; gdb -command $PWD/gdbCommands >& $procLog" >> $procCmdFile
-        echo "${node}$procCmdFile" >> $HOME/mpirun.schema
+        echo "${node}$procCmdFile" >> $PWD/mpirun.schema
     elif [ "$method" -eq 3 ]; then
         echo "$sourceFoam; cd $PWD; $exec $args >& $procLog" >> $procCmdFile
-        echo "${node}$procCmdFile" >> $HOME/mpirun.schema
+        echo "${node}$procCmdFile" >> $PWD/mpirun.schema
     elif [ "$method" -eq 4 ]; then
         echo "$sourceFoam; cd $PWD; $exec $args 2>&1 | tee $procLog; read dummy" >> $procCmdFile
-        echo "${node}xterm -font fixed -title 'processor'$proc $geom -e $procCmdFile" >> $HOME/mpirun.schema
+        echo "${node}xterm -font fixed -title 'processor'$proc $geom -e $procCmdFile" >> $PWD/mpirun.schema
     elif [ "$method" -eq 5 ]; then
         echo "$sourceFoam; cd $PWD; valgrind $exec $args; read dummy" >> $procCmdFile
-        echo "${node}xterm -font fixed -title 'processor'$proc $geom -e $procCmdFile" >> $HOME/mpirun.schema
+        echo "${node}xterm -font fixed -title 'processor'$proc $geom -e $procCmdFile" >> $PWD/mpirun.schema
     fi
 
     chmod +x $procCmdFile
@@ -176,10 +181,17 @@ do
     echo "    tail -f $procLog"
 done
 
-$ECHO "Constructed $HOME/mpirun.schema file. Press return to execute.\c"
+cmd=""
-read dummy
 if [ .$WM_MPLIB = .OPENMPI ]; then
-    mpirun -app $HOME/mpirun.schema </dev/null
+    cmd="mpirun -app $PWD/mpirun.schema </dev/null"
 elif [ .$WM_MPLIB = .LAM ]; then
-    mpirun $HOME/mpirun.schema </dev/null
+    cmd="mpirun $PWD/mpirun.schema </dev/null"
 fi
+
+echo "Constructed $PWD/mpirun.schema file."
+echo ""
+echo "    $cmd"
+echo ""
+$ECHO "Press return to execute.\c"
+read dummy
+exec $cmd

src/Allwmake

@@ -2,7 +2,8 @@
 cd ${0%/*} || exit 1    # run from this directory
 set -x
 
-( cd OpenFOAM && wmakeLnInclude . )
+wmakeLnInclude -f OpenFOAM
+wmakeLnInclude -f OSspecific/$WM_OS
 ( cd Pstream && ./Allwmake )
 
 wmake libo  OSspecific/$WM_OS

src/ODE/ODE/ODE.H

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

src/ODE/ODESolvers/KRR4/KRR4.C

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

src/ODE/ODESolvers/KRR4/KRR4.H

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

src/ODE/ODESolvers/SIBS/SIBS.H

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

src/ODE/ODESolvers/SIBS/SIMPR.C

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

src/ODE/ODESolvers/SIBS/polyExtrapolate.C

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

src/OSspecific/Unix/regularExpression.H

@@ -100,7 +100,7 @@ public:
     // Member functions
 
         //- Matches?
-        inline bool matches(const string& s)
+        inline bool matches(const string& s) const
         {
             size_t nmatch = 0;
             regmatch_t *pmatch = NULL;

src/OpenFOAM/Make/files

@@ -39,6 +39,7 @@ $(strings)/word/word.C
 $(strings)/word/wordIO.C
 $(strings)/fileName/fileName.C
 $(strings)/fileName/fileNameIO.C
+$(strings)/keyType/keyTypeIO.C
 
 primitives/random/Random.C
 
@@ -177,8 +178,9 @@ dimensionedTypes/dimensionedTensor/dimensionedTensor.C
 
 matrices/solution/solution.C
 
-scalarMatrix = matrices/scalarMatrix
+scalarMatrices = matrices/scalarMatrices
-$(scalarMatrix)/scalarMatrix.C
+$(scalarMatrices)/scalarMatrices.C
+$(scalarMatrices)/SVD/SVD.C
 
 LUscalarMatrix = matrices/LUscalarMatrix
 $(LUscalarMatrix)/LUscalarMatrix.C

src/OpenFOAM/db/IOstreams/IOstreams/Ostream.C

@@ -22,15 +22,31 @@ 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 "word.H"
 #include "Ostream.H"
 #include "token.H"
+#include "keyType.H"
+#include "IOstreams.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+//- Write keyType
+Foam::Ostream& Foam::Ostream::write(const keyType& s)
+{
+    // Write as word?
+    if (s.isWildCard())
+    {
+        return write(static_cast<const string&>(s));
+    }
+    else
+    {
+        return write(static_cast<const word&>(s));
+    }
+}
+
+
 //- Decrememt the indent level
 void Foam::Ostream::decrIndent()
 {
@@ -47,7 +63,7 @@ void Foam::Ostream::decrIndent()
 
 
 // Write the keyword to the Ostream followed by appropriate indentation
-Foam::Ostream& Foam::Ostream::writeKeyword(const Foam::word& keyword)
+Foam::Ostream& Foam::Ostream::writeKeyword(const Foam::keyType& keyword)
 {
     indent();
     write(keyword);

src/OpenFOAM/db/IOstreams/IOstreams/Ostream.H

@@ -35,6 +35,7 @@ Description
 #define Ostream_H
 
 #include "IOstream.H"
+#include "keyType.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -105,6 +106,9 @@ public:
             //- Write word
             virtual Ostream& write(const word&) = 0;
 
+            //- Write keyType
+            virtual Ostream& write(const keyType&);
+
             //- Write string
             virtual Ostream& write(const string&) = 0;
 
@@ -146,7 +150,7 @@ public:
 
             //- Write the keyword to the Ostream followed by
             //  appropriate indentation
-            Ostream& writeKeyword(const word& keyword);
+            Ostream& writeKeyword(const keyType& keyword);
 
 
         // Stream state functions

src/OpenFOAM/db/IOstreams/Pstreams/PstreamCommsStruct.C

@@ -22,8 +22,6 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
 #include "Pstream.H"

src/OpenFOAM/db/dictionary/dictionary.C

@@ -34,6 +34,72 @@ defineTypeNameAndDebug(Foam::dictionary, 0);
 
 const Foam::dictionary Foam::dictionary::null;
 
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+bool Foam::dictionary::findInWildcards
+(
+    const bool wildCardMatch,
+    const word& Keyword,
+    DLList<entry*>::const_iterator& wcLink,
+    DLList<autoPtr<regularExpression> >::const_iterator& reLink
+) const
+{
+    if (wildCardEntries_.size() > 0)
+    {
+        //wcLink = wildCardEntries_.begin();
+        //reLink = wildCardRegexps_.end();
+
+        while (wcLink != wildCardEntries_.end())
+        {
+            if (!wildCardMatch && wcLink()->keyword() == Keyword)
+            {
+                return true;
+            }
+            else if (wildCardMatch && reLink()->matches(Keyword))
+            {
+                return true;
+            }
+
+            ++reLink;
+            ++wcLink;
+        }
+    }
+
+    return false;
+}
+
+
+bool Foam::dictionary::findInWildcards
+(
+    const bool wildCardMatch,
+    const word& Keyword,
+    DLList<entry*>::iterator& wcLink,
+    DLList<autoPtr<regularExpression> >::iterator& reLink
+)
+{
+    if (wildCardEntries_.size() > 0)
+    {
+        while (wcLink != wildCardEntries_.end())
+        {
+            if (!wildCardMatch && wcLink()->keyword() == Keyword)
+            {
+                return true;
+            }
+            else if (wildCardMatch && reLink()->matches(Keyword))
+            {
+                return true;
+            }
+
+            ++reLink;
+            ++wcLink;
+        }
+    }
+
+    return false;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::dictionary::dictionary()
@@ -60,6 +126,18 @@ Foam::dictionary::dictionary
     )
     {
         hashedEntries_.insert(iter().keyword(), &iter());
+
+        if (iter().keyword().isWildCard())
+        {
+            wildCardEntries_.insert(&iter());
+            wildCardRegexps_.insert
+            (
+                autoPtr<regularExpression>
+                (
+                    new regularExpression(iter().keyword())
+                )
+            );
+        }
     }
 }
 
@@ -81,6 +159,18 @@ Foam::dictionary::dictionary
     )
     {
         hashedEntries_.insert(iter().keyword(), &iter());
+
+        if (iter().keyword().isWildCard())
+        {
+            wildCardEntries_.insert(&iter());
+            wildCardRegexps_.insert
+            (
+                autoPtr<regularExpression>
+                (
+                    new regularExpression(iter().keyword())
+                )
+            );
+        }
     }
 }
 
@@ -133,13 +223,29 @@ bool Foam::dictionary::found(const word& keyword, bool recursive) const
     {
         return true;
     }
-    else if (recursive && &parent_ != &dictionary::null)
-    {
-        return parent_.found(keyword, recursive);
-    }
     else
     {
-        return false;
+        if (wildCardEntries_.size() > 0)
+        {
+            DLList<entry*>::const_iterator wcLink = wildCardEntries_.begin();
+            DLList<autoPtr<regularExpression> >::const_iterator reLink =
+                wildCardRegexps_.begin();
+
+            // Find in wildcards using regular expressions only
+            if (findInWildcards(true, keyword, wcLink, reLink))
+            {
+                return true;
+            }
+        }
+
+        if (recursive && &parent_ != &dictionary::null)
+        {
+            return parent_.found(keyword, recursive);
+        }
+        else
+        {
+            return false;
+        }
     }
 }
 
@@ -147,16 +253,31 @@ bool Foam::dictionary::found(const word& keyword, bool recursive) const
 const Foam::entry* Foam::dictionary::lookupEntryPtr
 (
     const word& keyword,
-    bool recursive
+    bool recursive,
+    bool wildCardMatch
 ) const
 {
     HashTable<entry*>::const_iterator iter = hashedEntries_.find(keyword);
 
     if (iter == hashedEntries_.end())
     {
+        if (wildCardMatch && wildCardEntries_.size() > 0)
+        {
+            DLList<entry*>::const_iterator wcLink =
+                wildCardEntries_.begin();
+            DLList<autoPtr<regularExpression> >::const_iterator reLink =
+                wildCardRegexps_.begin();
+
+            // Find in wildcards using regular expressions only
+            if (findInWildcards(wildCardMatch, keyword, wcLink, reLink))
+            {
+                return wcLink();
+            }
+        }
+
         if (recursive && &parent_ != &dictionary::null)
         {
-            return parent_.lookupEntryPtr(keyword, recursive);
+            return parent_.lookupEntryPtr(keyword, recursive, wildCardMatch);
         }
         else
         {
@@ -171,19 +292,34 @@ const Foam::entry* Foam::dictionary::lookupEntryPtr
 Foam::entry* Foam::dictionary::lookupEntryPtr
 (
     const word& keyword,
-    bool recursive
+    bool recursive,
+    bool wildCardMatch
 )
 {
     HashTable<entry*>::iterator iter = hashedEntries_.find(keyword);
 
     if (iter == hashedEntries_.end())
     {
+        if (wildCardMatch && wildCardEntries_.size() > 0)
+        {
+            DLList<entry*>::iterator wcLink =
+                wildCardEntries_.begin();
+            DLList<autoPtr<regularExpression> >::iterator reLink =
+                wildCardRegexps_.begin();
+            // Find in wildcards using regular expressions only
+            if (findInWildcards(wildCardMatch, keyword, wcLink, reLink))
+            {
+                return wcLink();
+            }
+        }
+
         if (recursive && &parent_ != &dictionary::null)
         {
             return const_cast<dictionary&>(parent_).lookupEntryPtr
             (
                 keyword,
-                recursive
+                recursive,
+                wildCardMatch
             );
         }
         else
@@ -199,16 +335,17 @@ Foam::entry* Foam::dictionary::lookupEntryPtr
 const Foam::entry& Foam::dictionary::lookupEntry
 (
     const word& keyword,
-    bool recursive
+    bool recursive,
+    bool wildCardMatch
 ) const
 {
-    const entry* entryPtr = lookupEntryPtr(keyword, recursive);
+    const entry* entryPtr = lookupEntryPtr(keyword, recursive, wildCardMatch);
 
     if (entryPtr == NULL)
     {
         FatalIOErrorIn
         (
-            "dictionary::lookupEntry(const word& keyword) const",
+            "dictionary::lookupEntry(const word&, bool, bool) const",
             *this
         )   << "keyword " << keyword << " is undefined in dictionary "
             << name()
@@ -222,16 +359,18 @@ const Foam::entry& Foam::dictionary::lookupEntry
 Foam::ITstream& Foam::dictionary::lookup
 (
     const word& keyword,
-    bool recursive
+    bool recursive,
+    bool wildCardMatch
 ) const
 {
-    return lookupEntry(keyword, recursive).stream();
+    return lookupEntry(keyword, recursive, wildCardMatch).stream();
 }
 
 
 bool Foam::dictionary::isDict(const word& keyword) const
 {
-    const entry* entryPtr = lookupEntryPtr(keyword);
+    // Find non-recursive with wildcards
+    const entry* entryPtr = lookupEntryPtr(keyword, false, true);
 
     if (entryPtr)
     {
@@ -246,7 +385,7 @@ bool Foam::dictionary::isDict(const word& keyword) const
 
 const Foam::dictionary* Foam::dictionary::subDictPtr(const word& keyword) const
 {
-    const entry* entryPtr = lookupEntryPtr(keyword);
+    const entry* entryPtr = lookupEntryPtr(keyword, false, true);
 
     if (entryPtr)
     {
@@ -261,7 +400,8 @@ const Foam::dictionary* Foam::dictionary::subDictPtr(const word& keyword) const
 
 const Foam::dictionary& Foam::dictionary::subDict(const word& keyword) const
 {
-    const entry* entryPtr = lookupEntryPtr(keyword);
+    const entry* entryPtr = lookupEntryPtr(keyword, false, true);
+
     if (entryPtr == NULL)
     {
         FatalIOErrorIn
@@ -278,7 +418,8 @@ const Foam::dictionary& Foam::dictionary::subDict(const word& keyword) const
 
 Foam::dictionary& Foam::dictionary::subDict(const word& keyword)
 {
-    entry* entryPtr = lookupEntryPtr(keyword);
+    entry* entryPtr = lookupEntryPtr(keyword, false, true);
+
     if (entryPtr == NULL)
     {
         FatalIOErrorIn
@@ -314,7 +455,10 @@ Foam::wordList Foam::dictionary::toc() const
 
 bool Foam::dictionary::add(entry* entryPtr, bool mergeEntry)
 {
-    HashTable<entry*>::iterator iter = hashedEntries_.find(entryPtr->keyword());
+    HashTable<entry*>::iterator iter = hashedEntries_.find
+    (
+        entryPtr->keyword()
+    );
 
     if (mergeEntry && iter != hashedEntries_.end())
     {
@@ -336,6 +480,19 @@ bool Foam::dictionary::add(entry* entryPtr, bool mergeEntry)
             if (hashedEntries_.insert(entryPtr->keyword(), entryPtr))
             {
                 entryPtr->name() = name_ + "::" + entryPtr->keyword();
+
+                if (entryPtr->keyword().isWildCard())
+                {
+                    wildCardEntries_.insert(entryPtr);
+                    wildCardRegexps_.insert
+                    (
+                        autoPtr<regularExpression>
+                        (
+                            new regularExpression(entryPtr->keyword())
+                        )
+                    );
+                }
+
                 return true;
             }
             else
@@ -356,6 +513,18 @@ bool Foam::dictionary::add(entry* entryPtr, bool mergeEntry)
         entryPtr->name() = name_ + "::" + entryPtr->keyword();
         IDLList<entry>::append(entryPtr);
 
+        if (entryPtr->keyword().isWildCard())
+        {
+            wildCardEntries_.insert(entryPtr);
+            wildCardRegexps_.insert
+            (
+                autoPtr<regularExpression>
+                (
+                    new regularExpression(entryPtr->keyword())
+                )
+            );
+        }
+
         return true;
     }
     else
@@ -376,27 +545,37 @@ void Foam::dictionary::add(const entry& e, bool mergeEntry)
     add(e.clone(*this).ptr(), mergeEntry);
 }
 
-void Foam::dictionary::add(const word& k, const word& w, bool overwrite)
+void Foam::dictionary::add(const keyType& k, const word& w, bool overwrite)
 {
     add(new primitiveEntry(k, token(w)), overwrite);
 }
 
-void Foam::dictionary::add(const word& k, const Foam::string& s, bool overwrite)
+void Foam::dictionary::add
+(
+    const keyType& k,
+    const Foam::string& s,
+    bool overwrite
+)
 {
     add(new primitiveEntry(k, token(s)), overwrite);
 }
 
-void Foam::dictionary::add(const word& k, const label l, bool overwrite)
+void Foam::dictionary::add(const keyType& k, const label l, bool overwrite)
 {
     add(new primitiveEntry(k, token(l)), overwrite);
 }
 
-void Foam::dictionary::add(const word& k, const scalar s, bool overwrite)
+void Foam::dictionary::add(const keyType& k, const scalar s, bool overwrite)
 {
     add(new primitiveEntry(k, token(s)), overwrite);
 }
 
-void Foam::dictionary::add(const word& k, const dictionary& d, bool mergeEntry)
+void Foam::dictionary::add
+(
+    const keyType& k,
+    const dictionary& d,
+    bool mergeEntry
+)
 {
     add(new dictionaryEntry(k, *this, d), mergeEntry);
 }
@@ -404,7 +583,7 @@ void Foam::dictionary::add(const word& k, const dictionary& d, bool mergeEntry)
 
 void Foam::dictionary::set(entry* entryPtr)
 {
-    entry* existingPtr = lookupEntryPtr(entryPtr->keyword());
+    entry* existingPtr = lookupEntryPtr(entryPtr->keyword(), false, true);
 
     // clear dictionary so merge acts like overwrite
     if (existingPtr && existingPtr->isDict())
@@ -420,7 +599,7 @@ void Foam::dictionary::set(const entry& e)
     set(e.clone(*this).ptr());
 }
 
-void Foam::dictionary::set(const word& k, const dictionary& d)
+void Foam::dictionary::set(const keyType& k, const dictionary& d)
 {
     set(new dictionaryEntry(k, *this, d));
 }
@@ -432,6 +611,19 @@ bool Foam::dictionary::remove(const word& Keyword)
 
     if (iter != hashedEntries_.end())
     {
+        // Delete from wildcards first
+        DLList<entry*>::iterator wcLink =
+            wildCardEntries_.begin();
+        DLList<autoPtr<regularExpression> >::iterator reLink =
+            wildCardRegexps_.begin();
+
+        // Find in wildcards using exact match only
+        if (findInWildcards(false, Keyword, wcLink, reLink))
+        {
+            wildCardEntries_.remove(wcLink);
+            wildCardRegexps_.remove(reLink);
+        }
+
         IDLList<entry>::remove(iter());
         delete iter();
         hashedEntries_.erase(iter);
@@ -447,8 +639,8 @@ bool Foam::dictionary::remove(const word& Keyword)
 
 bool Foam::dictionary::changeKeyword
 (
-    const word& oldKeyword,
+    const keyType& oldKeyword,
-    const word& newKeyword,
+    const keyType& newKeyword,
     bool forceOverwrite
 )
 {
@@ -466,6 +658,18 @@ bool Foam::dictionary::changeKeyword
         return false;
     }
 
+    if (iter()->keyword().isWildCard())
+    {
+        FatalErrorIn
+        (
+            "dictionary::changeKeyword(const word&, const word&, bool)"
+        )   << "Old keyword "<< oldKeyword
+            << " is a wildcard."
+            << "Wildcard replacement not yet implemented."
+            << exit(FatalError);
+    }
+
+
     HashTable<entry*>::iterator iter2 = hashedEntries_.find(newKeyword);
 
     // newKeyword already exists
@@ -473,14 +677,33 @@ bool Foam::dictionary::changeKeyword
     {
         if (forceOverwrite)
         {
+            if (iter2()->keyword().isWildCard())
+            {
+                // Delete from wildcards first
+                DLList<entry*>::iterator wcLink =
+                    wildCardEntries_.begin();
+                DLList<autoPtr<regularExpression> >::iterator reLink =
+                    wildCardRegexps_.begin();
+
+                // Find in wildcards using exact match only
+                if (findInWildcards(false, iter2()->keyword(), wcLink, reLink))
+                {
+                    wildCardEntries_.remove(wcLink);
+                    wildCardRegexps_.remove(reLink);
+                }
+            }
+
             IDLList<entry>::replace(iter2(), iter());
             delete iter2();
             hashedEntries_.erase(iter2);
+            
         }
         else
         {
-            WarningIn("dictionary::changeKeyword(const word&, const word&)")
+            WarningIn
-                << "cannot rename keyword "<< oldKeyword
+            (
+                "dictionary::changeKeyword(const word&, const word&, bool)"
+            )   << "cannot rename keyword "<< oldKeyword
                 << " to existing keyword " << newKeyword
                 << " in dictionary " << name() << endl;
             return false;
@@ -493,6 +716,18 @@ bool Foam::dictionary::changeKeyword
     hashedEntries_.erase(oldKeyword);
     hashedEntries_.insert(newKeyword, iter());
 
+    if (newKeyword.isWildCard())
+    {
+        wildCardEntries_.insert(iter());
+        wildCardRegexps_.insert
+        (
+            autoPtr<regularExpression>
+            (
+                new regularExpression(newKeyword)
+            )
+        );
+    }
+
     return true;
 }
 
@@ -579,6 +814,7 @@ void Foam::dictionary::operator=(const dictionary& dict)
 
     // Create clones of the entries in the given dictionary
     // resetting the parentDict to this dictionary
+
     for
     (
         IDLList<entry>::const_iterator iter = dict.begin();
@@ -586,17 +822,7 @@ void Foam::dictionary::operator=(const dictionary& dict)
         ++iter
     )
     {
-        IDLList<entry>::append(iter().clone(*this).ptr());
+        add(iter().clone(*this).ptr());
-    }
-
-    for
-    (
-        IDLList<entry>::iterator iter = begin();
-        iter != end();
-        ++iter
-    )
-    {
-        hashedEntries_.insert(iter().keyword(), &iter());
     }
 }
 

src/OpenFOAM/db/dictionary/dictionary.H

@@ -27,7 +27,12 @@ Class
 
 Description
     A list of keyword definitions, which are a keyword followed by any number
-    of values (e.g. words and numbers).
+    of values (e.g. words and numbers). The keywords can represent wildcards
+    which are matched using Posix regular expressions. The general order for
+    searching is
+    - exact match
+    - wildcard match (in reverse order)
+    - optional recursion into subdictionaries
 
     The dictionary class is the base class for IOdictionary.
     It also serves as a bootstrap dictionary for the objectRegistry data
@@ -49,11 +54,13 @@ SourceFiles
 
 #include "entry.H"
 #include "IDLList.H"
+#include "DLList.H"
 #include "fileName.H"
 #include "ITstream.H"
 #include "HashTable.H"
 #include "wordList.H"
 #include "className.H"
+#include "regularExpression.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -80,12 +87,40 @@ class dictionary
         //- Dictionary name
         fileName name_;
 
-        //- HashTable of the enries held on the DL-list for quick lookup
+        //- HashTable of the entries held on the DL-list for quick lookup
         HashTable<entry*> hashedEntries_;
 
         //- Parent dictionary
         const dictionary& parent_;
 
+        //- Wildcard entries
+        DLList<entry*> wildCardEntries_;
+
+        //- Wildcard precompiled regex
+        DLList<autoPtr<regularExpression> > wildCardRegexps_;
+
+   // Private Member Functions
+
+        //- Search wildcard table either for exact match or for regular
+        //  expression match.
+        bool findInWildcards
+        (
+            const bool wildCardMatch,
+            const word& Keyword,
+            DLList<entry*>::const_iterator& wcLink,
+            DLList<autoPtr<regularExpression> >::const_iterator& reLink
+        ) const;
+
+        //- Search wildcard table either for exact match or for regular
+        //  expression match.
+        bool findInWildcards
+        (
+            const bool wildCardMatch,
+            const word& Keyword,
+            DLList<entry*>::iterator& wcLink,
+            DLList<autoPtr<regularExpression> >::iterator& reLink
+        );
+
 
 public:
 
@@ -181,24 +216,44 @@ public:
 
             //- Find and return an entry data stream pointer if present
             //  otherwise return NULL.
-            //  If recursive search parent dictionaries
+            //  If recursive search parent dictionaries. If wildCardMatch
+            //  use wildcards.
             const entry* lookupEntryPtr
             (
-                const word&, bool recursive=false
+                const word&,
+                bool recursive,
+                bool wildCardMatch
             ) const;
 
             //- Find and return an entry data stream pointer for manipulation
             //  if present otherwise return NULL.
-            //  If recursive search parent dictionaries
+            //  If recursive search parent dictionaries. If wildCardMatch
-            entry* lookupEntryPtr(const word&, bool recursive=false);
+            //  use wildcards.
+            entry* lookupEntryPtr
+            (
+                const word&,
+                bool recursive,
+                bool wildCardMatch
+            );
 
             //- Find and return an entry data stream if present otherwise error.
-            //  If recursive search parent dictionaries
+            //  If recursive search parent dictionaries. If wildCardMatch
-            const entry& lookupEntry(const word&, bool recursive=false) const;
+            //  use wildcards.
+            const entry& lookupEntry
+            (
+                const word&,
+                bool recursive,
+                bool wildCardMatch
+            ) const;
 
             //- Find and return an entry data stream
             //  If recursive search parent dictionaries
-            ITstream& lookup(const word&, bool recursive=false) const;
+            ITstream& lookup
+            (
+                const word&,
+                bool recursive=false,
+                bool wildCardMatch=true
+            ) const;
 
             //- Find and return a T,
             //  if not found return the given default value
@@ -208,7 +263,8 @@ public:
             (
                 const word&,
                 const T&,
-                bool recursive=false
+                bool recursive=false,
+                bool wildCardMatch=true
             ) const;
 
             //- Find and return a T, if not found return the given
@@ -219,7 +275,8 @@ public:
             (
                 const word&,
                 const T&,
-                bool recursive=false
+                bool recursive=false,
+                bool wildCardMatch=true
             );
 
             //- Find an entry if present, and assign to T
@@ -229,7 +286,8 @@ public:
             (
                 const word&,
                 T&,
-                bool recursive=false
+                bool recursive=false,
+                bool wildCardMatch=true
             ) const;
 
             //- Check if entry is a sub-dictionary
@@ -248,7 +306,6 @@ public:
             //- Return the table of contents
             wordList toc() const;
 
-
         // Editing
 
             //- Add a new entry
@@ -263,25 +320,25 @@ public:
 
             //- Add a word entry
             //  optionally overwrite an existing entry
-            void add(const word& keyword, const word&, bool overwrite=false);
+            void add(const keyType&, const word&, bool overwrite=false);
 
             //- Add a string entry
             //  optionally overwrite an existing entry
-            void add(const word& keyword, const string&, bool overwrite=false);
+            void add(const keyType&, const string&, bool overwrite=false);
 
             //- Add a label entry
             //  optionally overwrite an existing entry
-            void add(const word& keyword, const label, bool overwrite=false);
+            void add(const keyType&, const label, bool overwrite=false);
 
             //- Add a scalar entry
             //  optionally overwrite an existing entry
-            void add (const word& keyword, const scalar, bool overwrite=false);
+            void add (const keyType&, const scalar, bool overwrite=false);
 
             //- Add a dictionary entry
             //  optionally merge with an existing sub-dictionary
             void add
             (
-                const word& keyword,
+                const keyType& keyword,
                 const dictionary&,
                 bool mergeEntry=false
             );
@@ -289,7 +346,7 @@ public:
             //- Add a T entry
             //  optionally overwrite an existing entry
             template<class T>
-            void add(const word& keyword, const T&, bool overwrite=false);
+            void add(const keyType& keyword, const T&, bool overwrite=false);
 
             //- Assign a new entry, overwrite any existing entry
             void set(entry*);
@@ -298,11 +355,11 @@ public:
             void set(const entry&);
 
             //- Assign a dictionary entry, overwrite any existing entry
-            void set(const word& keyword, const dictionary&);
+            void set(const keyType& keyword, const dictionary&);
 
             //- Assign a T entry, overwrite any existing entry
             template<class T>
-            void set(const word& keyword, const T&);
+            void set(const keyType& keyword, const T&);
 
             //- Remove an entry specified by keyword
             bool remove(const word& keyword);
@@ -311,8 +368,8 @@ public:
             //  optionally forcing overwrite of an existing entry
             bool changeKeyword
             (
-                const word& oldKeyword,
+                const keyType& oldKeyword,
-                const word& newKeyword,
+                const keyType& newKeyword,
                 bool forceOverwrite = false
             );
 
@@ -361,11 +418,13 @@ public:
 
 // Global Operators
 
-//- Combine dictionaries starting from the entries in dict1 and then including those from dict2.
+//- Combine dictionaries starting from the entries in dict1 and then including
+//  those from dict2.
 //  Warn, but do not overwrite the entries from dict1.
 dictionary operator+(const dictionary& dict1, const dictionary& dict2);
 
-//- Combine dictionaries starting from the entries in dict1 and then including those from dict2.
+//- Combine dictionaries starting from the entries in dict1 and then including
+//  those from dict2.
 //  Do not overwrite the entries from dict1.
 dictionary operator|(const dictionary& dict1, const dictionary& dict2);
 

src/OpenFOAM/db/dictionary/dictionaryEntry/dictionaryEntry.C

@@ -30,7 +30,7 @@ License
 
 Foam::dictionaryEntry::dictionaryEntry
 (
-    const word& key,
+    const keyType& key,
     const dictionary& parentDict,
     const dictionary& dict
 )

src/OpenFOAM/db/dictionary/dictionaryEntry/dictionaryEntry.H

@@ -79,7 +79,7 @@ public:
         //- Construct from the keyword, parent dictionary and a Istream
         dictionaryEntry
         (
-            const word& keyword,
+            const keyType& keyword,
             const dictionary& parentDict,
             Istream& is
         );
@@ -87,7 +87,7 @@ public:
         //- Construct from the keyword, parent dictionary and a dictionary
         dictionaryEntry
         (
-            const word& keyword,
+            const keyType& keyword,
             const dictionary& parentDict,
             const dictionary& dict
         );

src/OpenFOAM/db/dictionary/dictionaryEntry/dictionaryEntryIO.C

@@ -27,7 +27,9 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
+#include "keyType.H"
 #include "dictionaryEntry.H"
+#include "IOstreams.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -43,14 +45,14 @@ Foam::dictionaryEntry::dictionaryEntry
     is.fatalCheck
     (
         "dictionaryEntry::dictionaryEntry"
-        "(Istream& is, const dictionary& parentDict)"
+        "(const dictionary& parentDict, Istream& is)"
     );
 }
 
 
 Foam::dictionaryEntry::dictionaryEntry
 (
-    const word& key,
+    const keyType& key,
     const dictionary& parentDict,
     Istream& is
 )
@@ -63,7 +65,7 @@ Foam::dictionaryEntry::dictionaryEntry
     is.fatalCheck
     (
         "dictionaryEntry::dictionaryEntry"
-        "(const word& keyword, const dictionary& parentDict, Istream& is)"
+        "(const keyType& keyword, const dictionary& parentDict, Istream& is)"
     );
 }
 

src/OpenFOAM/db/dictionary/dictionaryIO.C

@@ -71,7 +71,7 @@ bool Foam::dictionary::substituteKeyword(const word& keyword)
     word varName = keyword(1, keyword.size()-1);
 
     // lookup the variable name in the given dictionary....
-    const entry* ePtr = lookupEntryPtr(varName, true);
+    const entry* ePtr = lookupEntryPtr(varName, true, true);
 
     // ...if defined insert its entries into this dictionary...
     if (ePtr != NULL)
@@ -137,6 +137,8 @@ Foam::Istream& Foam::operator>>(Istream& is, dictionary& dict)
 
     dict.clear();
     dict.hashedEntries_.clear();
+    dict.wildCardEntries_.clear();
+    dict.wildCardRegexps_.clear();
     dict.read(is);
 
     return is;

src/OpenFOAM/db/dictionary/dictionaryTemplates.C

@@ -34,10 +34,11 @@ T Foam::dictionary::lookupOrDefault
 (
     const word& keyword,
     const T& deflt,
-    bool recursive
+    bool recursive,
+    bool wildCardMatch
 ) const
 {
-    const entry* entryPtr = lookupEntryPtr(keyword, recursive);
+    const entry* entryPtr = lookupEntryPtr(keyword, recursive, wildCardMatch);
 
     if (entryPtr == NULL)
     {
@@ -55,10 +56,11 @@ T Foam::dictionary::lookupOrAddDefault
 (
     const word& keyword,
     const T& deflt,
-    bool recursive
+    bool recursive,
+    bool wildCardMatch
 )
 {
-    const entry* entryPtr = lookupEntryPtr(keyword, recursive);
+    const entry* entryPtr = lookupEntryPtr(keyword, recursive, wildCardMatch);
 
     if (entryPtr == NULL)
     {
@@ -77,10 +79,11 @@ bool Foam::dictionary::readIfPresent
 (
     const word& k,
     T& val,
-    bool recursive
+    bool recursive,
+    bool wildCardMatch
 ) const
 {
-    const entry* entryPtr = lookupEntryPtr(k, recursive);
+    const entry* entryPtr = lookupEntryPtr(k, recursive, wildCardMatch);
 
     if (entryPtr == NULL)
     {
@@ -95,16 +98,17 @@ bool Foam::dictionary::readIfPresent
 
 
 template<class T>
-void Foam::dictionary::add(const word& k, const T& t, bool overwrite)
+void Foam::dictionary::add(const keyType& k, const T& t, bool overwrite)
 {
     add(new primitiveEntry(k, t), overwrite);
 }
 
 
 template<class T>
-void Foam::dictionary::set(const word& k, const T& t)
+void Foam::dictionary::set(const keyType& k, const T& t)
 {
     set(new primitiveEntry(k, t));
 }
 
+
 // ************************************************************************* //

src/OpenFOAM/db/dictionary/entry/entry.C

@@ -30,7 +30,7 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::entry::entry(const word& keyword)
+Foam::entry::entry(const keyType& keyword)
 :
     keyword_(keyword)
 {}

src/OpenFOAM/db/dictionary/entry/entry.H

@@ -42,6 +42,7 @@ SourceFiles
 #ifndef entry_H
 #define entry_H
 
+#include "keyType.H"
 #include "IDLList.H"
 #include "fileName.H"
 #include "autoPtr.H"
@@ -70,13 +71,13 @@ class entry
     // Private data
 
         //- Keyword of entry
-        word keyword_;
+        keyType keyword_;
 
 
     // Private Member Functions
 
         //- Get the next valid keyword otherwise return false
-        static bool getKeyword(word& keyword, Istream& is);
+        static bool getKeyword(keyType& keyword, Istream& is);
 
 
 public:
@@ -84,7 +85,7 @@ public:
     // Constructors
 
         //- Construct from keyword
-        entry(const word& keyword);
+        entry(const keyType& keyword);
 
         //- Construct as copy
         entry(const entry&);
@@ -116,13 +117,13 @@ public:
     // Member functions
 
         //- Return keyword
-        const word& keyword() const
+        const keyType& keyword() const
         {
             return keyword_;
         }
 
         //- Return non-const access to keyword
-        word& keyword()
+        keyType& keyword()
         {
             return keyword_;
         }

src/OpenFOAM/db/dictionary/entry/entryIO.C

@@ -32,7 +32,7 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-bool Foam::entry::getKeyword(word& keyword, Istream& is)
+bool Foam::entry::getKeyword(keyType& keyword, Istream& is)
 {
     token keywordToken;
 
@@ -57,6 +57,12 @@ bool Foam::entry::getKeyword(word& keyword, Istream& is)
         keyword = keywordToken.wordToken();
         return true;
     }
+    else if (keywordToken.isString())
+    {
+        // Enable wildcards
+        keyword = keywordToken.stringToken();
+        return true;
+    }
     // If it is the end of the dictionary or file return false...
     else if (keywordToken == token::END_BLOCK || is.eof())
     {
@@ -67,7 +73,7 @@ bool Foam::entry::getKeyword(word& keyword, Istream& is)
     {
         cerr<< "--> FOAM Warning : " << std::endl
             << "    From function "
-            << "entry::getKeyword(word& keyword, Istream& is)" << std::endl
+            << "entry::getKeyword(keyType& keyword, Istream& is)" << std::endl
             << "    in file " << __FILE__
             << " at line " << __LINE__ << std::endl
             << "    Reading " << is.name().c_str() << std::endl
@@ -84,7 +90,7 @@ bool Foam::entry::New(dictionary& parentDict, Istream& is)
 {
     is.fatalCheck("entry::New(const dictionary& parentDict, Istream& is)");
 
-    word keyword;
+    keyType keyword;
 
     // Get the next keyword and if invalid return false
     if (!getKeyword(keyword, is))
@@ -115,7 +121,13 @@ bool Foam::entry::New(dictionary& parentDict, Istream& is)
             // Deal with duplicate entries
             bool mergeEntry = false;
 
-            entry* existingPtr = parentDict.lookupEntryPtr(keyword);
+            // See (using exact match) if entry already present
+            entry* existingPtr = parentDict.lookupEntryPtr
+            (
+                keyword,
+                false,
+                false
+            );
             if (existingPtr)
             {
                 if (functionEntries::inputModeEntry::overwrite())
@@ -158,7 +170,7 @@ Foam::autoPtr<Foam::entry> Foam::entry::New(Istream& is)
 {
     is.fatalCheck("entry::New(Istream& is)");
 
-    word keyword;
+    keyType keyword;
 
     // Get the next keyword and if invalid return false
     if (!getKeyword(keyword, is))

src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntry.C

@@ -29,7 +29,7 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::primitiveEntry::primitiveEntry(const word& key, const ITstream& tokens)
+Foam::primitiveEntry::primitiveEntry(const keyType& key, const ITstream& tokens)
 :
     entry(key),
     ITstream(tokens)
@@ -38,7 +38,7 @@ Foam::primitiveEntry::primitiveEntry(const word& key, const ITstream& tokens)
 }
 
 
-Foam::primitiveEntry::primitiveEntry(const word& keyword, const token& t)
+Foam::primitiveEntry::primitiveEntry(const keyType& keyword, const token& t)
 :
     entry(keyword),
     ITstream(keyword, tokenList(1, t))
@@ -47,7 +47,7 @@ Foam::primitiveEntry::primitiveEntry(const word& keyword, const token& t)
 
 Foam::primitiveEntry::primitiveEntry
 (
-    const word& keyword,
+    const keyType& keyword,
     const tokenList& tokens
 )
 :

src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntry.H

@@ -108,23 +108,23 @@ public:
     // Constructors
 
         //- Construct from keyword and a Istream
-        primitiveEntry(const word& keyword, Istream&);
+        primitiveEntry(const keyType& keyword, Istream&);
 
         //- Construct from keyword, parent dictionary and a Istream
-        primitiveEntry(const word& keyword, const dictionary&, Istream&);
+        primitiveEntry(const keyType& keyword, const dictionary&, Istream&);
 
         //- Construct from keyword and a ITstream
-        primitiveEntry(const word& keyword, const ITstream&);
+        primitiveEntry(const keyType& keyword, const ITstream&);
 
         //- Construct from keyword and a token
-        primitiveEntry(const word&, const token&);
+        primitiveEntry(const keyType&, const token&);
 
         //- Construct from keyword and a tokenList
-        primitiveEntry(const word&, const tokenList&);
+        primitiveEntry(const keyType&, const tokenList&);
 
         //- Construct from keyword and a T
         template<class T>
-        primitiveEntry(const word&, const T&);
+        primitiveEntry(const keyType&, const T&);
 
         autoPtr<entry> clone(const dictionary&) const
         {

src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntryIO.C

@@ -81,7 +81,7 @@ bool Foam::primitiveEntry::expandVariable
     word varName = w(1, w.size()-1);
 
     // lookup the variable name in the given dictionary....
-    const entry* ePtr = dict.lookupEntryPtr(varName, true);
+    const entry* ePtr = dict.lookupEntryPtr(varName, true, true);
 
     // ...if defined insert its tokens into this 
     if (ePtr != NULL)
@@ -218,7 +218,7 @@ void Foam::primitiveEntry::readEntry(const dictionary& dict, Istream& is)
 
 Foam::primitiveEntry::primitiveEntry
 (
-    const word& key,
+    const keyType& key,
     const dictionary& dict,
     Istream& is
 )
@@ -236,7 +236,7 @@ Foam::primitiveEntry::primitiveEntry
 }
 
 
-Foam::primitiveEntry::primitiveEntry(const word& key, Istream& is)
+Foam::primitiveEntry::primitiveEntry(const keyType& key, Istream& is)
 :
     entry(key),
     ITstream

src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntryTemplates.C

@@ -30,7 +30,7 @@ License
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class T>
-Foam::primitiveEntry::primitiveEntry(const word& keyword, const T& t)
+Foam::primitiveEntry::primitiveEntry(const keyType& keyword, const T& t)
 :
     entry(keyword),
     ITstream(keyword, tokenList(10))

src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C

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

src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H

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

src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C

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

src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H

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

src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C

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

src/OpenFOAM/matrices/Matrix/Matrix.C

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

src/OpenFOAM/matrices/Matrix/Matrix.H

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

src/OpenFOAM/matrices/Matrix/MatrixI.H

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

src/OpenFOAM/matrices/Matrix/MatrixIO.C

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

src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H

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

src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H

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

src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H

@@ -0,0 +1,97 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 and rows
+        //  and value for all elements.
+        //  It checks that m == n.
+        inline SquareMatrix(const label m, const label n, const Type&);
+
+        //- Construct from Istream.
+        inline SquareMatrix(Istream&);
+
+        //- Clone
+        inline autoPtr<SquareMatrix<Type> > clone() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#   include "SquareMatrixI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H

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

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

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

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

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

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

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

src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C

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

src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H

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

src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C

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

src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.C

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

src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C

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

src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H

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

src/OpenFOAM/primitives/strings/keyType/keyType.H

@@ -23,88 +23,107 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::scalarMatrix
+    Foam::keyType
 
 Description
-    Foam::scalarMatrix
+    A class for handling keywords in dictionaries.
+
+    A keyType is the keyword of a dictionary. It differs from word in that
+    it accepts wildcards.
 
 SourceFiles
-    scalarMatrix.C
+    keyType.C
+    keyTypeIO.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef scalarMatrix_H
+#ifndef keyType_H
-#define scalarMatrix_H
+#define keyType_H
 
-#include "Matrix.H"
+#include "word.H"
-#include "scalarField.H"
-#include "labelList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
+// Forward declaration of classes
+class Istream;
+class Ostream;
+
+
 /*---------------------------------------------------------------------------*\
-                           Class scalarMatrix Declaration
+                           Class keyType Declaration
 \*---------------------------------------------------------------------------*/
 
-class scalarMatrix
+class keyType
 :
-    public Matrix<scalar>
+    public word
 {
+    // Private member data
+
+        bool isWildCard_;
+
+    // Private Member Functions
+
+        //- Disallow assignments where we cannot determine string/word type
+        void operator=(const std::string&);
 
 public:
 
+
     // Constructors
 
         //- Construct null
-        scalarMatrix();
+        inline keyType();
 
-        //- Construct given size
+        //- Construct as copy
-        scalarMatrix(const label);
+        inline keyType(const keyType& s);
 
-        //- Construct from Matrix<scalar>
+        //- Construct as copy of word
-        scalarMatrix(const Matrix<scalar>&);
+        inline keyType(const word& s);
+
+        //- Construct as copy of string. Expect it to be regular expression.
+        inline keyType(const string& s);
+
+        //- Construct as copy of character array
+        inline keyType(const char* s);
+
+        //- Construct as copy of std::string
+        inline keyType(const std::string& s, const bool isWildCard);
 
         //- Construct from Istream
-        scalarMatrix(Istream&);
+        keyType(Istream& is);
 
 
-    // Member Functions
+    // Member functions
 
-        //- Solve the matrix using Gaussian elimination with pivoting,
+        //- Is this character valid for a keyType
-        //  returning the solution in the source
+        inline static bool valid(char c);
-        template<class T>
-        static void solve(Matrix<scalar>& matrix, Field<T>& source);
 
-        //- Solve the matrix using Gaussian elimination with pivoting
+        //- Is the type a wildcard?
-        //  and return the solution
+        inline bool isWildCard() const;
-        template<class T>
-        void solve(Field<T>& psi, const Field<T>& source) const;
 
 
-        //- LU decompose the matrix with pivoting
+    // Member operators
-        static void LUDecompose
-        (
-            Matrix<scalar>& matrix,
-            labelList& pivotIndices
-        );
 
-        //- LU back-substitution with given source, returning the solution
+        // Assignment
-        //  in the source
-        template<class T>
-        static void LUBacksubstitute
-        (
-            const Matrix<scalar>& luMmatrix,
-            const labelList& pivotIndices,
-            Field<T>& source
-        );
 
-        //- Solve the matrix using LU decomposition with pivoting
+            inline void operator=(const keyType& s);
-        //  returning the LU form of the matrix and the solution in the source
+
-        template<class T>
+            //- Assign from regular expression.
-        static void LUsolve(Matrix<scalar>& matrix, Field<T>& source);
+            inline void operator=(const string& s);
+
+            inline void operator=(const word& s);
+
+            inline void operator=(const char*);
+
+
+    // IOstream operators
+
+        friend Istream& operator>>(Istream& is, keyType& w);
+
+        friend Ostream& operator<<(Ostream& os, const keyType& w);
 };
 
 
@@ -114,9 +133,7 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#ifdef NoRepository
+#include "keyTypeI.H"
-#   include "scalarMatrixTemplates.C"
-#endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/OpenFOAM/primitives/strings/keyType/keyTypeI.H

@@ -0,0 +1,132 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+//- Construct null
+inline Foam::keyType::keyType()
+:
+    word(),
+    isWildCard_(false)
+{}
+
+
+//- Construct as copy
+inline Foam::keyType::keyType(const keyType& s)
+:
+    word(s, false),
+    isWildCard_(s.isWildCard())
+{}
+
+
+//- Construct as copy of word
+inline Foam::keyType::keyType(const word& s)
+:
+    word(s, false),
+    isWildCard_(false)
+{}
+
+
+//- Construct as copy of string. Expect it to be regular expression
+inline Foam::keyType::keyType(const string& s)
+:
+    word(s, false),
+    isWildCard_(true)
+{}
+
+
+//- Construct as copy of character array
+inline Foam::keyType::keyType(const char* s)
+:
+    word(s, false),
+    isWildCard_(false)
+{}
+
+
+//- Construct as copy of std::string
+inline Foam::keyType::keyType
+(
+    const std::string& s,
+    const bool isWildCard
+)
+:
+    word(s, false),
+    isWildCard_(isWildCard)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline bool Foam::keyType::valid(char c)
+{
+    return c != '"';
+}
+
+
+bool Foam::keyType::isWildCard() const
+{
+    return isWildCard_;
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+inline void Foam::keyType::operator=(const keyType& s)
+{
+    // Bypass checking
+    string::operator=(s);
+    isWildCard_ = s.isWildCard();
+}
+
+
+inline void Foam::keyType::operator=(const word& s)
+{
+    word::operator=(s);
+    isWildCard_ = false;
+}
+
+
+inline void Foam::keyType::operator=(const string& s)
+{
+    // Bypass checking
+    string::operator=(s);
+    isWildCard_ = true;
+}
+
+
+inline void Foam::keyType::operator=(const char* s)
+{
+    // Bypass checking
+    string::operator=(s);
+    isWildCard_ = false;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/primitives/strings/keyType/keyTypeIO.C

@@ -0,0 +1,88 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Description
+    Istream constructor and IOstream operators for word.
+
+\*---------------------------------------------------------------------------*/
+
+#include "keyType.H"
+#include "IOstreams.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::keyType::keyType(Istream& is)
+:
+    word()
+{
+    is >> *this;
+}
+
+
+Foam::Istream& Foam::operator>>(Istream& is, keyType& w)
+{
+    token t(is);
+
+    if (!t.good())
+    {
+        is.setBad();
+        return is;
+    }
+
+    if (t.isWord())
+    {
+        w = t.wordToken();
+    }
+    else if (t.isString())
+    {
+        // Assign from string. Sets regular expression.
+        w = t.stringToken();
+    }
+    else
+    {
+        is.setBad();
+        FatalIOErrorIn("operator>>(Istream&, keyType&)", is)
+            << "wrong token type - expected word or string found "
+            << t.info()
+            << exit(FatalIOError);
+
+        return is;
+    }
+
+    // Check state of IOstream
+    is.check("Istream& operator>>(Istream&, keyType&)");
+
+    return is;
+}
+
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const keyType& w)
+{
+    os.write(w);
+    os.check("Ostream& operator<<(Ostream&, const keyType&)");
+    return os;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/primitives/strings/word/word.H

@@ -87,16 +87,21 @@ public:
         inline word(const word&);
 
         //- Construct as copy of character array
-        inline word(const char*);
+        inline word(const char*, const bool doStripInvalid = true);
 
         //- Construct as copy with a maximum number of characters
-        inline word(const char*, const size_type);
+        inline word
+        (
+            const char*,
+            const size_type,
+            const bool doStripInvalid
+        );
 
         //- Construct as copy of string
-        inline word(const string&);
+        inline word(const string&, const bool doStripInvalid = true);
 
         //- Construct as copy of std::string
-        inline word(const std::string&);
+        inline word(const std::string&, const bool doStripInvalid = true);
 
         //- Construct from Istream
         word(Istream&);

src/OpenFOAM/primitives/strings/word/wordI.H

@@ -65,34 +65,51 @@ inline Foam::word::word()
 {}
 
 
-inline Foam::word::word(const string& s)
+inline Foam::word::word(const string& s, const bool doStripInvalid)
 :
     string(s)
 {
-    stripInvalid();
+    if (doStripInvalid)
+    {
+        stripInvalid();
+    }
 }
 
 
-inline Foam::word::word(const std::string& s)
+inline Foam::word::word(const std::string& s, const bool doStripInvalid)
 :
     string(s)
 {
-    stripInvalid();
+    if (doStripInvalid)
+    {
+        stripInvalid();
+    }
 }
 
 
-inline Foam::word::word(const char* s)
+inline Foam::word::word(const char* s, const bool doStripInvalid)
 :
     string(s)
 {
-    stripInvalid();
+    if (doStripInvalid)
+    {
+        stripInvalid();
+    }
 }
 
-inline Foam::word::word(const char* s, const size_type n)
+inline Foam::word::word
+(
+    const char* s,
+    const size_type n,
+    const bool doStripInvalid
+)
 :
     string(s, n)
 {
-    stripInvalid();
+    if (doStripInvalid)
+    {
+        stripInvalid();
+    }
 }
 
 

src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C

@@ -826,6 +826,19 @@ void Foam::refinementSurfaces::findNearestIntersection
             }
         }
     }
+
+
+    // Make sure that if hit1 has hit something, hit2 will have at least the
+    // same point (due to tolerances it might miss its end point)
+    forAll(hit1, pointI)
+    {
+        if (hit1[pointI].hit() && !hit2[pointI].hit())
+        {
+            hit2[pointI] = hit1[pointI];
+            surface2[pointI] = surface1[pointI];
+            region2[pointI] = region1[pointI];
+        }
+    }
 }
 
 

src/finiteVolume/Make/files

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

src/finiteVolume/fvMesh/extendedStencil/extendedStencil.C

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

src/finiteVolume/fvMesh/extendedStencil/extendedStencil.H

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

src/finiteVolume/fvMesh/extendedStencil/extendedStencilTemplates.C

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

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

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

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

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

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

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

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

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

src/thermophysicalModels/barotropicCompressibilityModel/Chung/Chung.C

@@ -31,13 +31,16 @@ License
 
 namespace Foam
 {
-namespace compressibilityModels
+    namespace compressibilityModels
-{
+    {
-
+        defineTypeNameAndDebug(Chung, 0);
-defineTypeNameAndDebug(Chung, 0);
+        addToRunTimeSelectionTable
-addToRunTimeSelectionTable(barotropicCompressibilityModel, Chung, dictionary);
+        (
-
+            barotropicCompressibilityModel,
-}
+            Chung,
+            dictionary
+        );
+    }
 }
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -45,10 +48,11 @@ addToRunTimeSelectionTable(barotropicCompressibilityModel, Chung, dictionary);
 Foam::compressibilityModels::Chung::Chung
 (
     const dictionary& compressibilityProperties,
-    const volScalarField& gamma
+    const volScalarField& gamma,
+    const word& psiName
 )
 :
-    barotropicCompressibilityModel(compressibilityProperties, gamma),
+    barotropicCompressibilityModel(compressibilityProperties, gamma, psiName),
     psiv_(compressibilityProperties_.lookup("psiv")),
     psil_(compressibilityProperties_.lookup("psil")),
     rhovSat_(compressibilityProperties_.lookup("rhovSat")),

src/thermophysicalModels/barotropicCompressibilityModel/Chung/Chung.H

@@ -75,7 +75,8 @@ public:
         Chung
         (
             const dictionary& compressibilityProperties,
-            const volScalarField& gamma
+            const volScalarField& gamma,
+            const word& psiName = "psi"
         );
 
 

src/thermophysicalModels/barotropicCompressibilityModel/Wallis/Wallis.C

@@ -31,13 +31,16 @@ License
 
 namespace Foam
 {
-namespace compressibilityModels
+    namespace compressibilityModels
-{
+    {
-
+        defineTypeNameAndDebug(Wallis, 0);
-defineTypeNameAndDebug(Wallis, 0);
+        addToRunTimeSelectionTable
-addToRunTimeSelectionTable(barotropicCompressibilityModel, Wallis, dictionary);
+        (
-
+            barotropicCompressibilityModel,
-}
+            Wallis,
+            dictionary
+        );
+    }
 }
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -45,10 +48,11 @@ addToRunTimeSelectionTable(barotropicCompressibilityModel, Wallis, dictionary);
 Foam::compressibilityModels::Wallis::Wallis
 (
     const dictionary& compressibilityProperties,
-    const volScalarField& gamma
+    const volScalarField& gamma,
+    const word& psiName
 )
 :
-    barotropicCompressibilityModel(compressibilityProperties, gamma),
+    barotropicCompressibilityModel(compressibilityProperties, gamma, psiName),
     psiv_(compressibilityProperties_.lookup("psiv")),
     psil_(compressibilityProperties_.lookup("psil")),
     rhovSat_(compressibilityProperties_.lookup("rhovSat")),
@@ -62,8 +66,9 @@ Foam::compressibilityModels::Wallis::Wallis
 
 void Foam::compressibilityModels::Wallis::correct()
 {
-    psi_ = (gamma_*rhovSat_ + (scalar(1) - gamma_)*rholSat_)
+    psi_ =
-         *(gamma_*psiv_/rhovSat_ + (scalar(1) - gamma_)*psil_/rholSat_);
+        (gamma_*rhovSat_ + (scalar(1) - gamma_)*rholSat_)
+       *(gamma_*psiv_/rhovSat_ + (scalar(1) - gamma_)*psil_/rholSat_);
 }
 
 

src/thermophysicalModels/barotropicCompressibilityModel/Wallis/Wallis.H

@@ -75,7 +75,8 @@ public:
         Wallis
         (
             const dictionary& compressibilityProperties,
-            const volScalarField& gamma
+            const volScalarField& gamma,
+            const word& psiName = "psi"
         );
 
 

src/thermophysicalModels/barotropicCompressibilityModel/barotropicCompressibilityModel/barotropicCompressibilityModel.C

@@ -42,7 +42,8 @@ namespace Foam
 Foam::barotropicCompressibilityModel::barotropicCompressibilityModel
 (
     const dictionary& compressibilityProperties,
-    const volScalarField& gamma
+    const volScalarField& gamma,
+    const word& psiName
 )
 :
     compressibilityProperties_(compressibilityProperties),
@@ -50,12 +51,12 @@ Foam::barotropicCompressibilityModel::barotropicCompressibilityModel
     (
         IOobject
         (
-            "psi",
+            psiName,
             gamma.mesh().time().timeName(),
             gamma.mesh()
         ),
         gamma.mesh(),
-        dimensionedScalar("psi", dimensionSet(0, -2, 2, 0, 0), 0)
+        dimensionedScalar(psiName, dimensionSet(0, -2, 2, 0, 0), 0)
     ),
     gamma_(gamma)
 {}

src/thermophysicalModels/barotropicCompressibilityModel/barotropicCompressibilityModel/barotropicCompressibilityModel.H

@@ -97,9 +97,10 @@ public:
             dictionary,
             (
                 const dictionary& compressibilityProperties,
-                const volScalarField& gamma
+                const volScalarField& gamma,
+                const word& psiName
             ),
-            (compressibilityProperties, gamma)
+            (compressibilityProperties, gamma, psiName)
         );
 
 
@@ -109,7 +110,8 @@ public:
         static autoPtr<barotropicCompressibilityModel> New
         (
             const dictionary& compressibilityProperties,
-            const volScalarField& gamma
+            const volScalarField& gamma,
+            const word& psiName = "psi"
         );
 
 
@@ -119,7 +121,8 @@ public:
         barotropicCompressibilityModel
         (
             const dictionary& compressibilityProperties,
-            const volScalarField& gamma
+            const volScalarField& gamma,
+            const word& psiName = "psi"
         );
 
 

src/thermophysicalModels/barotropicCompressibilityModel/barotropicCompressibilityModel/newBarotropicCompressibilityModel.C

@@ -32,7 +32,8 @@ Foam::autoPtr<Foam::barotropicCompressibilityModel>
 Foam::barotropicCompressibilityModel::New
 (
     const dictionary& compressibilityProperties,
-    const volScalarField& gamma
+    const volScalarField& gamma,
+    const word& psiName
 )
 {
     word bcModelTypeName
@@ -60,7 +61,7 @@ Foam::barotropicCompressibilityModel::New
 
     return autoPtr<barotropicCompressibilityModel>
     (
-        cstrIter()(compressibilityProperties, gamma)
+        cstrIter()(compressibilityProperties, gamma, psiName)
     );
 }