Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

README.html

@@ -102,7 +102,7 @@ lang="en" xml:lang="en">
 <h2 id="sec-1"><span class="section-number-2">1</span> About OpenFOAM </h2>
 <div class="outline-text-2" id="text-1">
 
-<p>  OpenFOAM is a free, open source computational fluid dynamcis (CFD) software
+<p>  OpenFOAM is a free, open source computational fluid dynamics (CFD) software
   package released by the OpenFOAM Foundation. It has a large user base across
   most areas of engineering and science, from both commercial and academic
   organisations. OpenFOAM has an extensive range of features to solve anything

README.org

@@ -8,7 +8,7 @@
 # Copyright (c) 2011 OpenFOAM Foundation.
 
 * About OpenFOAM
-  OpenFOAM is a free, open source computational fluid dynamcis (CFD) software
+  OpenFOAM is a free, open source computational fluid dynamics (CFD) software
   package released by the OpenFOAM Foundation. It has a large user base across
   most areas of engineering and science, from both commercial and academic
   organisations. OpenFOAM has an extensive range of features to solve anything

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/EEqns.H

@@ -27,7 +27,7 @@
       - (fvc::ddt(alpha1) + fvc::div(alphaPhi1))*K1
 
       + (
-            he1.name() == "e"
+            he1.name() == thermo1.phasePropertyName("e")
           ? fvc::div(alphaPhi1, p)
           : -dalpha1pdt
         )/rho1
@@ -49,7 +49,7 @@
       - (fvc::ddt(alpha2) + fvc::div(alphaPhi2))*K2
 
       + (
-            he2.name() == "e"
+            he2.name() == thermo2.phasePropertyName("e")
           ? fvc::div(alphaPhi2, p)
           : -dalpha2pdt
         )/rho2

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H

@@ -1,3 +1,7 @@
+mrfZones.correctBoundaryVelocity(U1);
+mrfZones.correctBoundaryVelocity(U2);
+mrfZones.correctBoundaryVelocity(U);
+
 fvVectorMatrix U1Eqn(U1, U1.dimensions()*dimVol/dimTime);
 fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
 

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H

@@ -104,6 +104,35 @@ surfaceScalarField alphaPhi2("alphaPhi" + phase2Name, phi2);
             {
                 alphaPhi1 = alphaPhic1;
             }
+
+            /*
+            // Legacy semi-implicit and potentially unbounded form
+            fvScalarMatrix alpha1Eqn
+            (
+                fvm::ddt(alpha1)
+              + fvm::div(phic, alpha1, alphaScheme)
+              + fvm::div
+                (
+                    -fvc::flux(-phir, alpha2, alpharScheme),
+                    alpha1,
+                    alpharScheme
+                )
+             ==
+               fvm::Sp(Sp, alpha1) + Su
+            );
+
+            alpha1Eqn.relax();
+            alpha1Eqn.solve();
+
+            if (nAlphaSubCycles > 1)
+            {
+                alphaPhi1 += (runTime.deltaT()/totalDeltaT)*alpha1Eqn.flux();
+            }
+            else
+            {
+                alphaPhi1 = alpha1Eqn.flux();
+            }
+            */
         }
 
         if (g0.value() > 0.0)

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.C

@@ -26,6 +26,7 @@ License
 #include "phaseModel.H"
 #include "diameterModel.H"
 #include "fixedValueFvPatchFields.H"
+#include "slipFvPatchFields.H"
 #include "surfaceInterpolate.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -118,7 +119,11 @@ Foam::phaseModel::phaseModel
 
         forAll(U_.boundaryField(), i)
         {
-            if (isA<fixedValueFvPatchVectorField>(U_.boundaryField()[i]))
+            if
+            (
+                isA<fixedValueFvPatchVectorField>(U_.boundaryField()[i])
+             || isA<slipFvPatchVectorField>(U_.boundaryField()[i])
+            )
             {
                 phiTypes[i] = fixedValueFvPatchScalarField::typeName;
             }

applications/solvers/multiphase/multiphaseEulerFoam/UEqns.H

@@ -1,3 +1,5 @@
+#include "mrfZonesCorrectBCs.H"
+
 PtrList<fvVectorMatrix> UEqns(fluid.phases().size());
 autoPtr<multiphaseSystem::dragCoeffFields> dragCoeffs(fluid.dragCoeffs());
 

applications/solvers/multiphase/multiphaseEulerFoam/mrfZonesCorrectBCs.H

@@ -0,0 +1,6 @@
+    forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+    {
+        mrfZones.correctBoundaryVelocity(iter().U());
+    }
+
+    mrfZones.correctBoundaryVelocity(U);

applications/solvers/multiphase/multiphaseEulerFoam/phaseModel/phaseModel/phaseModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,6 +26,7 @@ License
 #include "phaseModel.H"
 #include "diameterModel.H"
 #include "fixedValueFvPatchFields.H"
+#include "slipFvPatchFields.H"
 #include "surfaceInterpolate.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -152,7 +153,11 @@ Foam::phaseModel::phaseModel
 
         forAll(U_.boundaryField(), i)
         {
-            if (isA<fixedValueFvPatchVectorField>(U_.boundaryField()[i]))
+            if
+            (
+                isA<fixedValueFvPatchVectorField>(U_.boundaryField()[i])
+             || isA<slipFvPatchVectorField>(U_.boundaryField()[i])
+            )
             {
                 phiTypes[i] = fixedValueFvPatchScalarField::typeName;
             }

applications/solvers/multiphase/twoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,6 +25,7 @@ License
 
 #include "phaseModel.H"
 #include "fixedValueFvPatchFields.H"
+#include "slipFvPatchFields.H"
 #include "surfaceInterpolate.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -114,7 +115,11 @@ Foam::phaseModel::phaseModel
 
         forAll(U_.boundaryField(), i)
         {
-            if (isA<fixedValueFvPatchVectorField>(U_.boundaryField()[i]))
+            if
+            (
+                isA<fixedValueFvPatchVectorField>(U_.boundaryField()[i])
+             || isA<slipFvPatchVectorField>(U_.boundaryField()[i])
+            )
             {
                 phiTypes[i] = fixedValueFvPatchScalarField::typeName;
             }

applications/test/Matrix/Test-Matrix.C

@@ -23,7 +23,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "SquareMatrix.H"
+#include "scalarMatrices.H"
 #include "vector.H"
 
 using namespace Foam;
@@ -49,15 +49,15 @@ int main(int argc, char *argv[])
     Info<< max(hmm) << endl;
     Info<< min(hmm) << endl;
 
-    SquareMatrix<scalar> hmm2(3, 1.0);
+    SquareMatrix<scalar> hmm2(3, 3, 1.0);
 
     hmm = hmm2;
 
     Info<< hmm << endl;
 
-    SquareMatrix<scalar> hmm3(Sin);
+    //SquareMatrix<scalar> hmm3(Sin);
 
-    Info<< hmm3 << endl;
+    //Info<< hmm3 << endl;
 
     SquareMatrix<scalar> hmm4;
 
@@ -70,7 +70,65 @@ int main(int argc, char *argv[])
     hmm4 = hmm5;
     Info<< hmm5 << endl;
 
-    Info<< "End\n" << endl;
+    {
+        scalarSymmetricSquareMatrix symmMatrix(3, 3, 0);
+
+        symmMatrix(0, 0) = 4;
+        symmMatrix(1, 0) = 12;
+        symmMatrix(1, 1) = 37;
+        symmMatrix(2, 0) = -16;
+        symmMatrix(2, 1) = -43;
+        symmMatrix(2, 2) = 98;
+
+        Info<< "Symmetric Square Matrix = " << symmMatrix << endl;
+
+        Info<< "Inverse = " << inv(symmMatrix) << endl;
+        Info<< "Determinant = " << det(symmMatrix) << endl;
+
+        scalarSymmetricSquareMatrix symmMatrix2(symmMatrix);
+        LUDecompose(symmMatrix2);
+
+        Info<< "Inverse = " << invDecomposed(symmMatrix2) << endl;
+        Info<< "Determinant = " << detDecomposed(symmMatrix2) << endl;
+
+        scalarDiagonalMatrix rhs(3, 0);
+        rhs[0] = 1;
+        rhs[1] = 2;
+        rhs[2] = 3;
+
+        LUsolve(symmMatrix, rhs);
+
+        Info<< "Decomposition = " << symmMatrix << endl;
+        Info<< "Solution = " << rhs << endl;
+    }
+
+    {
+        scalarSquareMatrix squareMatrix(3, 3, 0);
+
+        squareMatrix[0][0] = 4;
+        squareMatrix[0][1] = 12;
+        squareMatrix[0][2] = -16;
+        squareMatrix[1][0] = 12;
+        squareMatrix[1][1] = 37;
+        squareMatrix[1][2] = -43;
+        squareMatrix[2][0] = -16;
+        squareMatrix[2][1] = -43;
+        squareMatrix[2][2] = 98;
+
+        Info<< nl << "Square Matrix = " << squareMatrix << endl;
+
+        scalarDiagonalMatrix rhs(3, 0);
+        rhs[0] = 1;
+        rhs[1] = 2;
+        rhs[2] = 3;
+
+        LUsolve(squareMatrix, rhs);
+
+        Info<< "Decomposition = " << squareMatrix << endl;
+        Info<< "Solution = " << rhs << endl;
+    }
+
+    Info<< "\nEnd\n" << endl;
 
     return 0;
 }

etc/config/settings.csh

@@ -226,7 +226,7 @@ case ThirdParty:
         breaksw
     case Gcc47:
     case Gcc47++0x:
-        set gcc_version=gcc-4.7.0
+        set gcc_version=gcc-4.7.2
         set gmp_version=gmp-5.0.4
         set mpfr_version=mpfr-3.1.0
         set mpc_version=mpc-0.9

etc/config/settings.sh

@@ -247,7 +247,7 @@ OpenFOAM | ThirdParty)
         mpc_version=mpc-0.9
         ;;
     Gcc47 | Gcc47++0x)
-        gcc_version=gcc-4.7.0
+        gcc_version=gcc-4.7.2
         gmp_version=gmp-5.0.4
         mpfr_version=mpfr-3.1.0
         mpc_version=mpc-0.9

src/OpenFOAM/db/objectRegistry/objectRegistry.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -159,6 +159,10 @@ public:
             template<class Type>
             HashTable<const Type*> lookupClass(const bool strict = false) const;
 
+            //- Lookup and return all objects of the given Type
+            template<class Type>
+            HashTable<Type*> lookupClass(const bool strict = false);
+
             //- Is the named Type found?
             template<class Type>
             bool foundObject(const word& name) const;

src/OpenFOAM/db/objectRegistry/objectRegistryTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -76,6 +76,34 @@ Foam::HashTable<const Type*> Foam::objectRegistry::lookupClass
 }
 
 
+template<class Type>
+Foam::HashTable<Type*> Foam::objectRegistry::lookupClass
+(
+    const bool strict
+)
+{
+    HashTable<Type*> objectsOfClass(size());
+
+    forAllIter(HashTable<regIOobject*>, *this, iter)
+    {
+        if
+        (
+            (strict && isType<Type>(*iter()))
+         || (!strict && isA<Type>(*iter()))
+        )
+        {
+            objectsOfClass.insert
+            (
+                iter()->name(),
+                dynamic_cast<Type*>(iter())
+            );
+        }
+    }
+
+    return objectsOfClass;
+}
+
+
 template<class Type>
 bool Foam::objectRegistry::foundObject(const word& name) const
 {

src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -76,7 +76,7 @@ public:
 
     // Member functions
 
-        //- Invert the diaganol matrix and return itself
+        //- Invert the diagonal matrix and return itself
         DiagonalMatrix<Type>& invert();
 };
 

src/OpenFOAM/matrices/Matrix/Matrix.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -221,7 +221,6 @@ template<class Form, class Type> Form operator*
     const Matrix<Form, Type>&
 );
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.C

@@ -0,0 +1,98 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "SymmetricSquareMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::SymmetricSquareMatrix<Type> Foam::invDecomposed
+(
+    const SymmetricSquareMatrix<Type>& matrix
+)
+{
+    SymmetricSquareMatrix<Type> inv(matrix.n(), matrix.n(), 0.0);
+
+    for (label i = 0; i < matrix.n(); ++i)
+    {
+        inv[i][i] = 1.0/matrix[i][i];
+
+        for (label j = 0; j < i; ++j)
+        {
+            scalar sum = 0.0;
+
+            for (label k = j; k < i; k++)
+            {
+                sum -= matrix[i][k]*inv[k][j];
+            }
+
+            inv[i][j] = sum/matrix[i][i];
+        }
+    }
+
+    return inv.T()*inv;
+}
+
+
+template<class Type>
+Foam::SymmetricSquareMatrix<Type> Foam::inv
+(
+    const SymmetricSquareMatrix<Type>& matrix
+)
+{
+    SymmetricSquareMatrix<Type> matrixTmp(matrix);
+
+    LUDecompose(matrixTmp);
+
+    return invDecomposed(matrixTmp);
+}
+
+
+template<class Type>
+Foam::scalar Foam::detDecomposed(const SymmetricSquareMatrix<Type>& matrix)
+{
+    scalar diagProduct = 1.0;
+
+    for (label i = 0; i < matrix.n(); ++i)
+    {
+        diagProduct *= matrix[i][i];
+    }
+
+    return sqr(diagProduct);
+}
+
+
+template<class Type>
+Foam::scalar Foam::det(const SymmetricSquareMatrix<Type>& matrix)
+{
+    SymmetricSquareMatrix<Type> matrixTmp = matrix;
+
+    LUDecompose(matrixTmp);
+
+    return detDecomposed(matrixTmp);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H

@@ -0,0 +1,126 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::SymmetricSquareMatrix
+
+Description
+    A templated 2D square symmetric matrix of objects of \<T\>, where the
+    n x n matrix dimension is known and used for subscript bounds checking, etc.
+
+SourceFiles
+    SymmetricSquareMatrixI.H
+    SymmetricSquareMatrix.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef SymmetricSquareMatrix_H
+#define SymmetricSquareMatrix_H
+
+#include "SquareMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class SymmetricSquareMatrix Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class SymmetricSquareMatrix
+:
+    public Matrix<SymmetricSquareMatrix<Type>, Type>
+{
+
+public:
+
+    // Constructors
+
+        //- Null constructor.
+        inline SymmetricSquareMatrix();
+
+        //- Construct given number of rows/columns.
+        inline SymmetricSquareMatrix(const label n);
+
+        //- Construct with given number of rows/columns
+        inline SymmetricSquareMatrix(const label m, const label n);
+
+        //- Construct with given number of rows/columns
+        //  and value for all elements.
+        inline SymmetricSquareMatrix(const label m, const label n, const Type&);
+
+        //- Construct from Istream.
+        inline SymmetricSquareMatrix(Istream&);
+
+        //- Clone
+        inline autoPtr<SymmetricSquareMatrix<Type> > clone() const;
+
+
+        //- Return subscript-checked row of Matrix.
+        inline Type& operator()(const label r, const label c);
+
+        //- Return subscript-checked row of constant Matrix.
+        inline const Type& operator()(const label r, const label c) const;
+};
+
+
+// Global functions
+
+//- Return the LU decomposed SymmetricSquareMatrix inverse
+template<class Type>
+SymmetricSquareMatrix<Type> invDecomposed(const SymmetricSquareMatrix<Type>&);
+
+//- Return the SymmetricSquareMatrix inverse
+template<class Type>
+SymmetricSquareMatrix<Type> inv(const SymmetricSquareMatrix<Type>&);
+
+//- Return the LU decomposed SymmetricSquareMatrix det
+template<class Type>
+scalar detDecomposed(const SymmetricSquareMatrix<Type>&);
+
+//- Return the SymmetricSquareMatrix det
+template<class Type>
+scalar det(const SymmetricSquareMatrix<Type>&);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#   include "SymmetricSquareMatrixI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "SymmetricSquareMatrix.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrixI.H

@@ -0,0 +1,139 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix()
+:
+    Matrix<SymmetricSquareMatrix<Type>, Type>()
+{}
+
+
+template<class Type>
+inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix(const label n)
+:
+    Matrix<SymmetricSquareMatrix<Type>, Type>(n, n)
+{}
+
+
+template<class Type>
+inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix
+(
+    const label m,
+    const label n
+)
+:
+    Matrix<SymmetricSquareMatrix<Type>, Type>(m, n)
+{
+    if (m != n)
+    {
+        FatalErrorIn
+        (
+            "SymmetricSquareMatrix<Type>::SymmetricSquareMatrix"
+            "(const label m, const label n)"
+        )   << "m != n for constructing a symmetric square matrix"
+            << exit(FatalError);
+    }
+}
+
+
+template<class Type>
+inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix
+(
+    const label m,
+    const label n,
+    const Type& t
+)
+:
+    Matrix<SymmetricSquareMatrix<Type>, Type>(m, n, t)
+{
+    if (m != n)
+    {
+        FatalErrorIn
+        (
+            "SymmetricSquareMatrix<Type>::SymmetricSquareMatrix"
+            "(const label m, const label n, const Type&)"
+        )   << "m != n for constructing a symmetric square matrix"
+            << exit(FatalError);
+    }
+}
+
+
+template<class Type>
+inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix(Istream& is)
+:
+    Matrix<SymmetricSquareMatrix<Type>, Type>(is)
+{}
+
+
+template<class Type>
+inline Foam::autoPtr<Foam::SymmetricSquareMatrix<Type> >
+Foam::SymmetricSquareMatrix<Type>::clone() const
+{
+    return autoPtr<SymmetricSquareMatrix<Type> >
+    (
+        new SymmetricSquareMatrix<Type>(*this)
+    );
+}
+
+
+template<class Type>
+inline Type& Foam::SymmetricSquareMatrix<Type>::operator()
+(
+    const label r,
+    const label c
+)
+{
+    if (r > c)
+    {
+        return this->operator[](r)[c];
+    }
+    else
+    {
+        return this->operator[](c)[r];
+    }
+}
+
+
+template<class Type>
+inline const Type& Foam::SymmetricSquareMatrix<Type>::operator()
+(
+    const label r,
+    const label c
+) const
+{
+    if (r > c)
+    {
+        return this->operator[](r)[c];
+    }
+    else
+    {
+        return this->operator[](c)[r];
+    }
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -70,7 +70,7 @@ Foam::GAMGAgglomeration::GAMGAgglomeration
     const dictionary& controlDict
 )
 :
-    MeshObject<lduMesh, GAMGAgglomeration>(mesh),
+    MeshObject<lduMesh, Foam::GeometricMeshObject, GAMGAgglomeration>(mesh),
 
     maxLevels_(50),
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -58,7 +58,7 @@ class lduMatrix;
 
 class GAMGAgglomeration
 :
-    public MeshObject<lduMesh, GAMGAgglomeration>
+    public MeshObject<lduMesh, GeometricMeshObject, GAMGAgglomeration>
 {
 protected:
 

src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -134,6 +134,56 @@ void Foam::LUDecompose
 }
 
 
+void Foam::LUDecompose(scalarSymmetricSquareMatrix& matrix)
+{
+    // Store result in upper triangular part of matrix
+    label size = matrix.n();
+
+    // Set upper triangular parts to zero.
+    for (label j = 0; j < size; j++)
+    {
+        for (label k = j + 1; k < size; k++)
+        {
+            matrix[j][k] = 0.0;
+        }
+    }
+
+    for (label j = 0; j < size; j++)
+    {
+        scalar d = 0.0;
+
+        for (label k = 0; k < j; k++)
+        {
+            scalar s = 0.0;
+
+            for (label i = 0; i < k; i++)
+            {
+                s += matrix[i][k]*matrix[i][j];
+            }
+
+            s = (matrix[j][k] - s)/matrix[k][k];
+
+            matrix[k][j] = s;
+            matrix[j][k] = s;
+
+            d += sqr(s);
+        }
+
+        d = matrix[j][j] - d;
+
+        if (d < 0.0)
+        {
+            FatalErrorIn("Foam::LUDecompose(scalarSymmetricSquareMatrix&)")
+                << "Matrix is not symmetric positive-definite. Unable to "
+                << "decompose."
+                << abort(FatalError);
+        }
+
+        matrix[j][j] = sqrt(d);
+    }
+}
+
+
 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 
 void Foam::multiply

src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,6 +27,10 @@ Class
 Description
     Scalar matrices
 
+    LUDecompose for scalarSymmetricSquareMatrix implements the Cholesky
+    decomposition method from JAMA, a public-domain library developed at NIST,
+    available at http://math.nist.gov/tnt/index.html
+
 SourceFiles
     scalarMatrices.C
     scalarMatricesTemplates.C
@@ -38,6 +42,7 @@ SourceFiles
 
 #include "RectangularMatrix.H"
 #include "SquareMatrix.H"
+#include "SymmetricSquareMatrix.H"
 #include "DiagonalMatrix.H"
 #include "scalarField.H"
 #include "labelList.H"
@@ -49,21 +54,22 @@ namespace Foam
 
 typedef RectangularMatrix<scalar> scalarRectangularMatrix;
 typedef SquareMatrix<scalar> scalarSquareMatrix;
+typedef SymmetricSquareMatrix<scalar> scalarSymmetricSquareMatrix;
 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);
+void solve(scalarSquareMatrix& matrix, List<Type>& source);
 
 //- Solve the matrix using Gaussian elimination with pivoting
 //  and return the solution
 template<class Type>
 void solve
 (
-    Field<Type>& psi,
+    List<Type>& psi,
     const scalarSquareMatrix& matrix,
-    const Field<Type>& source
+    const List<Type>& source
 );
 
 //- LU decompose the matrix with pivoting
@@ -73,6 +79,9 @@ void LUDecompose
     labelList& pivotIndices
 );
 
+//- LU decompose the matrix into a lower (L) and upper (U) part. U = L.T()
+void LUDecompose(scalarSymmetricSquareMatrix& matrix);
+
 //- LU back-substitution with given source, returning the solution
 //  in the source
 template<class Type>
@@ -80,13 +89,28 @@ void LUBacksubstitute
 (
     const scalarSquareMatrix& luMmatrix,
     const labelList& pivotIndices,
-    Field<Type>& source
+    List<Type>& source
+);
+
+//- LU back-substitution with given source, returning the solution
+//  in the source. Specialised for symmetric square matrices that have been
+//  decomposed into LU where U = L.T() as pivoting is not required
+template<class Type>
+void LUBacksubstitute
+(
+    const scalarSymmetricSquareMatrix& luMmatrix,
+    List<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 LUsolve(scalarSquareMatrix& matrix, List<Type>& source);
+
+//- Solve the matrix using LU decomposition returning the LU form of the matrix
+//  and the solution in the source, where U = L.T()
+template<class Type>
+void LUsolve(scalarSymmetricSquareMatrix& matrix, List<Type>& source);
 
 template<class Form, class Type>
 void multiply

src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,6 +25,7 @@ License
 
 #include "scalarMatrices.H"
 #include "Swap.H"
+#include "ListOps.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -32,7 +33,7 @@ template<class Type>
 void Foam::solve
 (
     scalarSquareMatrix& tmpMatrix,
-    Field<Type>& sourceSol
+    List<Type>& sourceSol
 )
 {
     label n = tmpMatrix.n();
@@ -103,9 +104,9 @@ void Foam::solve
 template<class Type>
 void Foam::solve
 (
-    Field<Type>& psi,
+    List<Type>& psi,
     const scalarSquareMatrix& matrix,
-    const Field<Type>& source
+    const List<Type>& source
 )
 {
     scalarSquareMatrix tmpMatrix = matrix;
@@ -119,7 +120,7 @@ void Foam::LUBacksubstitute
 (
     const scalarSquareMatrix& luMatrix,
     const labelList& pivotIndices,
-    Field<Type>& sourceSol
+    List<Type>& sourceSol
 )
 {
     label n = luMatrix.n();
@@ -163,11 +164,57 @@ void Foam::LUBacksubstitute
 }
 
 
+template<class Type>
+void Foam::LUBacksubstitute
+(
+    const scalarSymmetricSquareMatrix& luMatrix,
+    List<Type>& sourceSol
+)
+{
+    label n = luMatrix.n();
+
+    label ii = 0;
+
+    for (register label i=0; i<n; i++)
+    {
+        Type sum = sourceSol[i];
+        const scalar* __restrict__ luMatrixi = luMatrix[i];
+
+        if (ii != 0)
+        {
+            for (label j=ii-1; j<i; j++)
+            {
+                sum -= luMatrixi[j]*sourceSol[j];
+            }
+        }
+        else if (sum != pTraits<Type>::zero)
+        {
+            ii = i+1;
+        }
+
+        sourceSol[i] = sum/luMatrixi[i];
+    }
+
+    for (register label i=n-1; i>=0; 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];
+        }
+
+        sourceSol[i] = sum/luMatrixi[i];
+    }
+}
+
+
 template<class Type>
 void Foam::LUsolve
 (
     scalarSquareMatrix& matrix,
-    Field<Type>& sourceSol
+    List<Type>& sourceSol
 )
 {
     labelList pivotIndices(matrix.n());
@@ -176,6 +223,18 @@ void Foam::LUsolve
 }
 
 
+template<class Type>
+void Foam::LUsolve
+(
+    scalarSymmetricSquareMatrix& matrix,
+    List<Type>& sourceSol
+)
+{
+    LUDecompose(matrix);
+    LUBacksubstitute(matrix, sourceSol);
+}
+
+
 template<class Form, class Type>
 void Foam::multiply
 (

src/OpenFOAM/meshes/MeshObject/MeshObject.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,26 +28,18 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-template<class Mesh, class Type>
+template<class Mesh, template<class> class MeshObjectType, class Type>
-Foam::MeshObject<Mesh, Type>::MeshObject(const Mesh& mesh)
+Foam::MeshObject<Mesh, MeshObjectType, Type>::MeshObject(const Mesh& mesh)
 :
-    regIOobject
+    MeshObjectType<Mesh>(Type::typeName, mesh.thisDb()),
-    (
-        IOobject
-        (
-            Type::typeName,
-            mesh.thisDb().instance(),
-            mesh.thisDb()
-        )
-    ),
     mesh_(mesh)
 {}
 
 
 // * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
 
-template<class Mesh, class Type>
+template<class Mesh, template<class> class MeshObjectType, class Type>
-const Type& Foam::MeshObject<Mesh, Type>::New
+const Type& Foam::MeshObject<Mesh, MeshObjectType, Type>::New
 (
     const Mesh& mesh
 )
@@ -67,14 +59,14 @@ const Type& Foam::MeshObject<Mesh, Type>::New
     }
     else
     {
-        return store(new Type(mesh));
+        return regIOobject::store(new Type(mesh));
     }
 }
 
 
-template<class Mesh, class Type>
+template<class Mesh, template<class> class MeshObjectType, class Type>
 template<class Data1>
-const Type& Foam::MeshObject<Mesh, Type>::New
+const Type& Foam::MeshObject<Mesh, MeshObjectType, Type>::New
 (
     const Mesh& mesh,
     const Data1& d
@@ -95,14 +87,14 @@ const Type& Foam::MeshObject<Mesh, Type>::New
     }
     else
     {
-        return store(new Type(mesh, d));
+        return regIOobject::store(new Type(mesh, d));
     }
 }
 
 
-template<class Mesh, class Type>
+template<class Mesh, template<class> class MeshObjectType, class Type>
 template<class Data1, class Data2>
-const Type& Foam::MeshObject<Mesh, Type>::New
+const Type& Foam::MeshObject<Mesh, MeshObjectType, Type>::New
 (
     const Mesh& mesh,
     const Data1& d1,
@@ -124,14 +116,14 @@ const Type& Foam::MeshObject<Mesh, Type>::New
     }
     else
     {
-        return store(new Type(mesh, d1, d2));
+        return regIOobject::store(new Type(mesh, d1, d2));
     }
 }
 
 
-template<class Mesh, class Type>
+template<class Mesh, template<class> class MeshObjectType, class Type>
 template<class Data1, class Data2, class Data3>
-const Type& Foam::MeshObject<Mesh, Type>::New
+const Type& Foam::MeshObject<Mesh, MeshObjectType, Type>::New
 (
     const Mesh& mesh,
     const Data1& d1,
@@ -154,14 +146,14 @@ const Type& Foam::MeshObject<Mesh, Type>::New
     }
     else
     {
-        return store(new Type(mesh, d1, d2, d3));
+        return regIOobject::store(new Type(mesh, d1, d2, d3));
     }
 }
 
 
-template<class Mesh, class Type>
+template<class Mesh, template<class> class MeshObjectType, class Type>
 template<class Data1, class Data2, class Data3, class Data4>
-const Type& Foam::MeshObject<Mesh, Type>::New
+const Type& Foam::MeshObject<Mesh, MeshObjectType, Type>::New
 (
     const Mesh& mesh,
     const Data1& d1,
@@ -185,15 +177,15 @@ const Type& Foam::MeshObject<Mesh, Type>::New
     }
     else
     {
-        return store(new Type(mesh, d1, d2, d3, d4));
+        return regIOobject::store(new Type(mesh, d1, d2, d3, d4));
     }
 }
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * //
 
-template<class Mesh, class Type>
+template<class Mesh, template<class> class MeshObjectType, class Type>
-bool Foam::MeshObject<Mesh, Type>::Delete(const Mesh& mesh)
+bool Foam::MeshObject<Mesh, MeshObjectType, Type>::Delete(const Mesh& mesh)
 {
     if
     (
@@ -221,10 +213,79 @@ bool Foam::MeshObject<Mesh, Type>::Delete(const Mesh& mesh)
 }
 
 
-template<class Mesh, class Type>
+template<class Mesh, template<class> class MeshObjectType, class Type>
-Foam::MeshObject<Mesh, Type>::~MeshObject()
+Foam::MeshObject<Mesh, MeshObjectType, Type>::~MeshObject()
 {
-    release();
+    MeshObjectType<Mesh>::release();
+}
+
+
+template<class Mesh>
+void Foam::meshObject::movePoints(objectRegistry& obr)
+{
+    HashTable<GeometricMeshObject<Mesh>*> meshObjects
+    (
+        obr.lookupClass<GeometricMeshObject<Mesh> >()
+    );
+
+    forAllIter
+    (
+        typename HashTable<GeometricMeshObject<Mesh>*>,
+        meshObjects,
+        iter
+    )
+    {
+        if (isA<MoveableMeshObject<Mesh> >(*iter()))
+        {
+            dynamic_cast<MoveableMeshObject<Mesh>*>(iter())->movePoints();
+        }
+        else
+        {
+            obr.checkOut(*iter());
+        }
+    }
+}
+
+
+template<class Mesh>
+void Foam::meshObject::updateMesh(objectRegistry& obr, const mapPolyMesh& mpm)
+{
+    HashTable<GeometricMeshObject<Mesh>*> meshObjects
+    (
+        obr.lookupClass<GeometricMeshObject<Mesh> >()
+    );
+
+    forAllIter
+    (
+        typename HashTable<GeometricMeshObject<Mesh>*>,
+        meshObjects,
+        iter
+    )
+    {
+        if (isA<UpdateableMeshObject<Mesh> >(*iter()))
+        {
+            dynamic_cast<UpdateableMeshObject<Mesh>*>(iter())->updateMesh(mpm);
+        }
+        else
+        {
+            obr.checkOut(*iter());
+        }
+    }
+}
+
+
+template<class Mesh, template<class> class MeshObjectType>
+void Foam::meshObject::clear(objectRegistry& obr)
+{
+    HashTable<MeshObjectType<Mesh>*> meshObjects
+    (
+        obr.lookupClass<MeshObjectType<Mesh> >()
+    );
+
+    forAllIter(typename HashTable<MeshObjectType<Mesh>*>, meshObjects, iter)
+    {
+        obr.checkOut(*iter());
+    }
 }
 
 

src/OpenFOAM/meshes/MeshObject/MeshObject.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,9 +25,36 @@ Class
     Foam::MeshObject
 
 Description
-    Templated abstract base-class for dynamic mesh objects used to automate
+    Templated abstract base-class for optional mesh objects used to automate
     their allocation to the mesh database and the mesh-modifier event-loop.
 
+    MeshObject is templated on the type of mesh it is allocated to, the type of
+    the mesh object (TopologicalMeshObject, GeometricMeshObject,
+    MoveableMeshObject, UpdateableMeshObject) and the type of the actual object
+    it is created for example:
+
+    class leastSquaresVectors
+    :
+        public MeshObject<fvMesh, MoveableMeshObject, leastSquaresVectors>
+    {
+    .
+    .
+    .
+        //- Delete the least square vectors when the mesh moves
+        virtual bool movePoints();
+    };
+
+    MeshObject types:
+
+    TopologicalMeshObject: mesh object to be deleted on topology change
+    GeometricMeshObject: mesh object to be deleted on geometry change
+    MoveableMeshObject: mesh object to be updated in movePoints
+    UpdateableMeshObject: mesh object to be updated in updateMesh or movePoints
+
+    Note that movePoints must be provided for MeshObjects of type
+    MoveableMeshObject and both movePoints and updateMesh functions must exist
+    provided for MeshObjects of type UpdateableMeshObject.
+
 SourceFiles
     MeshObject.C
 
@@ -37,21 +64,24 @@ SourceFiles
 #define MeshObject_H
 
 #include "regIOobject.H"
+#include "objectRegistry.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
+// Forward declarations
+class mapPolyMesh;
+
 /*---------------------------------------------------------------------------*\
                          Class MeshObject Declaration
 \*---------------------------------------------------------------------------*/
 
-
+template<class Mesh, template<class> class MeshObjectType, class Type>
-template<class Mesh, class Type>
 class MeshObject
 :
-    public regIOobject
+    public MeshObjectType<Mesh>
 {
 
 protected:
@@ -124,6 +154,121 @@ public:
 };
 
 
+/*---------------------------------------------------------------------------*\
+                           Class meshObject Declaration
+\*---------------------------------------------------------------------------*/
+
+class meshObject
+:
+    public regIOobject
+{
+public:
+
+    // Constructors
+
+        meshObject(const word& typeName, const objectRegistry& obr)
+        :
+            regIOobject
+            (
+                IOobject
+                (
+                    typeName,
+                    obr.instance(),
+                    obr
+                )
+            )
+        {}
+
+
+    // Static member functions
+
+        template<class Mesh>
+        static void movePoints(objectRegistry&);
+
+        template<class Mesh>
+        static void updateMesh(objectRegistry&, const mapPolyMesh&);
+
+        template<class Mesh, template<class> class MeshObjectType>
+        static void clear(objectRegistry&);
+};
+
+
+/*---------------------------------------------------------------------------*\
+                    Class TopologicalMeshObject Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Mesh>
+class TopologicalMeshObject
+:
+    public meshObject
+{
+public:
+
+    TopologicalMeshObject(const word& typeName, const objectRegistry& obr)
+    :
+        meshObject(typeName, obr)
+    {}
+};
+
+
+/*---------------------------------------------------------------------------*\
+                    Class GeometricMeshObject Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Mesh>
+class GeometricMeshObject
+:
+    public TopologicalMeshObject<Mesh>
+{
+public:
+
+    GeometricMeshObject(const word& typeName, const objectRegistry& obr)
+    :
+        TopologicalMeshObject<Mesh>(typeName, obr)
+    {}
+};
+
+
+/*---------------------------------------------------------------------------*\
+                    Class MoveableMeshObject Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Mesh>
+class MoveableMeshObject
+:
+    public GeometricMeshObject<Mesh>
+{
+public:
+
+    MoveableMeshObject(const word& typeName, const objectRegistry& obr)
+    :
+        GeometricMeshObject<Mesh>(typeName, obr)
+    {}
+
+    virtual bool movePoints() = 0;
+};
+
+
+/*---------------------------------------------------------------------------*\
+                    Class UpdateableMeshObject Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Mesh>
+class UpdateableMeshObject
+:
+    public MoveableMeshObject<Mesh>
+{
+public:
+
+    UpdateableMeshObject(const word& typeName, const objectRegistry& obr)
+    :
+        MoveableMeshObject<Mesh>(typeName, obr)
+    {}
+
+    virtual void updateMesh(const mapPolyMesh& mpm) = 0;
+};
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/meshes/pointMesh/pointMesh.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -72,7 +72,7 @@ void Foam::pointMesh::mapFields(const mapPolyMesh& mpm)
 
 Foam::pointMesh::pointMesh(const polyMesh& pMesh)
 :
-    MeshObject<polyMesh, pointMesh>(pMesh),
+    MeshObject<polyMesh, Foam::UpdateableMeshObject, pointMesh>(pMesh),
     GeoMesh<polyMesh>(pMesh),
     boundary_(*this, pMesh.boundaryMesh())
 {
@@ -88,7 +88,7 @@ Foam::pointMesh::pointMesh(const polyMesh& pMesh)
 }
 
 
-void Foam::pointMesh::movePoints(const pointField& newPoints)
+bool Foam::pointMesh::movePoints()
 {
     if (debug)
     {
@@ -96,7 +96,9 @@ void Foam::pointMesh::movePoints(const pointField& newPoints)
             << "Moving points." << endl;
     }
 
-    boundary_.movePoints(newPoints);
+    boundary_.movePoints(GeoMesh<polyMesh>::mesh_.points());
+
+    return true;
 }
 
 

src/OpenFOAM/meshes/pointMesh/pointMesh.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -48,7 +48,7 @@ namespace Foam
 
 class pointMesh
 :
-    public MeshObject<polyMesh, pointMesh>,
+    public MeshObject<polyMesh, UpdateableMeshObject, pointMesh>,
     public GeoMesh<polyMesh>
 {
     // Permanent data
@@ -121,7 +121,7 @@ public:
         // Mesh motion
 
             //- Move points, returns volumes swept by faces in motion
-            void movePoints(const pointField&);
+            bool movePoints();
 
             //- Update the mesh corresponding to given map
             void updateMesh(const mapPolyMesh& mpm);

src/OpenFOAM/meshes/polyMesh/polyMesh.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,30 +26,24 @@ License
 #include "polyMesh.H"
 #include "Time.H"
 #include "cellIOList.H"
-#include "SubList.H"
 #include "wedgePolyPatch.H"
 #include "emptyPolyPatch.H"
 #include "globalMeshData.H"
 #include "processorPolyPatch.H"
-#include "OSspecific.H"
 #include "polyMeshTetDecomposition.H"
 #include "indexedOctree.H"
 #include "treeDataCell.H"
-#include "SubField.H"
+#include "MeshObject.H"
 
-#include "pointMesh.H"
-#include "Istream.H"
-#include "Ostream.H"
-#include "simpleRegIOobject.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-defineTypeNameAndDebug(polyMesh, 0);
+    defineTypeNameAndDebug(polyMesh, 0);
 
-word polyMesh::defaultRegion = "region0";
+    word polyMesh::defaultRegion = "region0";
-word polyMesh::meshSubDir = "polyMesh";
+    word polyMesh::meshSubDir = "polyMesh";
 }
 
 
@@ -1162,21 +1156,7 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
     geometricD_ = Vector<label>::zero;
     solutionD_ = Vector<label>::zero;
 
-
+    meshObject::movePoints<polyMesh>(*this);
-    // Hack until proper callbacks. Below are all the polyMeh MeshObjects with a
-    // movePoints function.
-
-    // pointMesh
-    if (thisDb().foundObject<pointMesh>(pointMesh::typeName))
-    {
-        const_cast<pointMesh&>
-        (
-            thisDb().lookupObject<pointMesh>
-            (
-                pointMesh::typeName
-            )
-        ).movePoints(points_);
-    }
 
     const_cast<Time&>(time()).functionObjects().movePoints(*this);
 

src/OpenFOAM/meshes/polyMesh/polyMeshClear.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,8 +26,7 @@ License
 #include "polyMesh.H"
 #include "primitiveMesh.H"
 #include "globalMeshData.H"
-#include "pointMesh.H"
+#include "MeshObject.H"
-#include "Time.H"
 #include "indexedOctree.H"
 #include "treeDataCell.H"
 
@@ -61,6 +60,8 @@ void Foam::polyMesh::clearGeom()
             << endl;
     }
 
+    meshObject::clear<polyMesh, GeometricMeshObject>(*this);
+
     primitiveMesh::clearGeom();
 
     boundary_.clearGeom();
@@ -101,6 +102,8 @@ void Foam::polyMesh::clearAddressing()
             << endl;
     }
 
+    meshObject::clear<polyMesh, TopologicalMeshObject>(*this);
+
     primitiveMesh::clearAddressing();
 
     // parallelData depends on the processorPatch ordering so force
@@ -118,6 +121,7 @@ void Foam::polyMesh::clearAddressing()
 
     // Remove the stored tet base points
     tetBasePtIsPtr_.clear();
+
     // Remove the cell tree
     cellTreePtr_.clear();
 }
@@ -132,8 +136,6 @@ void Foam::polyMesh::clearPrimitives()
     owner_.setSize(0);
     neighbour_.setSize(0);
 
-    pointMesh::Delete(*this);
-
     clearedPrimitives_ = true;
 }
 
@@ -142,8 +144,6 @@ void Foam::polyMesh::clearOut()
 {
     clearGeom();
     clearAddressing();
-
-    pointMesh::Delete(*this);
 }
 
 

src/OpenFOAM/meshes/polyMesh/polyMeshUpdate.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -91,25 +91,12 @@ void Foam::polyMesh::updateMesh(const mapPolyMesh& mpm)
         }
     }
 
+    meshObject::updateMesh<polyMesh>(*this, mpm);
+
     // Reset valid directions (could change by faces put into empty patches)
     geometricD_ = Vector<label>::zero;
     solutionD_ = Vector<label>::zero;
 
-
-    // Hack until proper callbacks. Below are all the polyMesh-MeshObjects.
-
-    // pointMesh
-    if (thisDb().foundObject<pointMesh>(pointMesh::typeName))
-    {
-        const_cast<pointMesh&>
-        (
-            thisDb().lookupObject<pointMesh>
-            (
-                pointMesh::typeName
-            )
-        ).updateMesh(mpm);
-    }
-
     const_cast<Time&>(time()).functionObjects().updateMesh(mpm);
 }
 

src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -302,11 +302,11 @@ Foam::dynamicRefineFvMesh::refine
             Pout<< "Found " << masterFaces.size() << " split faces " << endl;
         }
 
-        HashTable<const surfaceScalarField*> fluxes
+        HashTable<surfaceScalarField*> fluxes
         (
             lookupClass<surfaceScalarField>()
         );
-        forAllConstIter(HashTable<const surfaceScalarField*>, fluxes, iter)
+        forAllIter(HashTable<surfaceScalarField*>, fluxes, iter)
         {
             if (!correctFluxes_.found(iter.key()))
             {
@@ -336,7 +336,7 @@ Foam::dynamicRefineFvMesh::refine
                     << endl;
             }
 
-            surfaceScalarField& phi = const_cast<surfaceScalarField&>(*iter());
+            surfaceScalarField& phi = *iter();
             const surfaceScalarField phiU
             (
                 fvc::interpolate
@@ -519,11 +519,11 @@ Foam::dynamicRefineFvMesh::unrefine
         const labelList& reversePointMap = map().reversePointMap();
         const labelList& reverseFaceMap = map().reverseFaceMap();
 
-        HashTable<const surfaceScalarField*> fluxes
+        HashTable<surfaceScalarField*> fluxes
         (
             lookupClass<surfaceScalarField>()
         );
-        forAllConstIter(HashTable<const surfaceScalarField*>, fluxes, iter)
+        forAllIter(HashTable<surfaceScalarField*>, fluxes, iter)
         {
             if (!correctFluxes_.found(iter.key()))
             {
@@ -553,7 +553,7 @@ Foam::dynamicRefineFvMesh::unrefine
                     << endl;
             }
 
-            surfaceScalarField& phi = const_cast<surfaceScalarField&>(*iter());
+            surfaceScalarField& phi = *iter();
             surfaceScalarField::GeometricBoundaryField& bphi =
                 phi.boundaryField();
 

src/dynamicMesh/fvMeshDistribute/fvMeshDistributeTemplates.C

@@ -66,7 +66,7 @@ void Foam::fvMeshDistribute::saveBoundaryFields
 
     HashTable<const fldType*> flds
     (
-        mesh_.objectRegistry::lookupClass<fldType>()
+        static_cast<const fvMesh&>(mesh_).objectRegistry::lookupClass<fldType>()
     );
 
     bflds.setSize(flds.size());
@@ -97,7 +97,7 @@ void Foam::fvMeshDistribute::mapBoundaryFields
 
     typedef GeometricField<T, fvsPatchField, Mesh> fldType;
 
-    HashTable<const fldType*> flds
+    HashTable<fldType*> flds
     (
         mesh_.objectRegistry::lookupClass<fldType>()
     );
@@ -110,15 +110,11 @@ void Foam::fvMeshDistribute::mapBoundaryFields
 
     label fieldI = 0;
 
-    forAllConstIter(typename HashTable<const fldType*>, flds, iter)
+    forAllIter(typename HashTable<fldType*>, flds, iter)
     {
-        const fldType& fld = *iter();
+        fldType& fld = *iter();
         typename fldType::GeometricBoundaryField& bfld =
-            const_cast<typename fldType::GeometricBoundaryField&>
+            fld.boundaryField();
-            (
-                fld.boundaryField()
-            );
-
 
         const FieldField<fvsPatchField, T>& oldBfld = oldBflds[fieldI++];
 
@@ -156,20 +152,17 @@ void Foam::fvMeshDistribute::initPatchFields
     const typename GeoField::value_type& initVal
 )
 {
-    HashTable<const GeoField*> flds
+    HashTable<GeoField*> flds
     (
         mesh_.objectRegistry::lookupClass<GeoField>()
     );
 
-    forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
+    forAllIter(typename HashTable<GeoField*>, flds, iter)
     {
-        const GeoField& fld = *iter();
+        GeoField& fld = *iter();
 
         typename GeoField::GeometricBoundaryField& bfld =
-            const_cast<typename GeoField::GeometricBoundaryField&>
+            fld.boundaryField();
-            (
-                fld.boundaryField()
-            );
 
         forAll(bfld, patchI)
         {
@@ -186,16 +179,15 @@ void Foam::fvMeshDistribute::initPatchFields
 template<class GeoField>
 void Foam::fvMeshDistribute::correctBoundaryConditions()
 {
-    HashTable<const GeoField*> flds
+    HashTable<GeoField*> flds
     (
         mesh_.objectRegistry::lookupClass<GeoField>()
     );
 
-    forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
+    forAllIter(typename HashTable<GeoField*>, flds, iter)
     {
         const GeoField& fld = *iter();
-
+        fld.correctBoundaryConditions();
-        const_cast<GeoField&>(fld).correctBoundaryConditions();
     }
 }
 

src/dynamicMesh/fvMeshTools/fvMeshToolsTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,20 +38,17 @@ void Foam::fvMeshTools::addPatchFields
     const typename GeoField::value_type& defaultPatchValue
 )
 {
-    HashTable<const GeoField*> flds
+    HashTable<GeoField*> flds
     (
         mesh.objectRegistry::lookupClass<GeoField>()
     );
 
-    forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
+    forAllIter(typename HashTable<GeoField*>, flds, iter)
     {
-        const GeoField& fld = *iter();
+        GeoField& fld = *iter();
 
         typename GeoField::GeometricBoundaryField& bfld =
-            const_cast<typename GeoField::GeometricBoundaryField&>
+            fld.boundaryField();
-            (
-                fld.boundaryField()
-            );
 
         label sz = bfld.size();
         bfld.setSize(sz+1);
@@ -95,20 +92,17 @@ void Foam::fvMeshTools::setPatchFields
     const dictionary& patchFieldDict
 )
 {
-    HashTable<const GeoField*> flds
+    HashTable<GeoField*> flds
     (
         mesh.objectRegistry::lookupClass<GeoField>()
     );
 
-    forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
+    forAllIter(typename HashTable<GeoField*>, flds, iter)
     {
-        const GeoField& fld = *iter();
+        GeoField& fld = *iter();
 
         typename GeoField::GeometricBoundaryField& bfld =
-            const_cast<typename GeoField::GeometricBoundaryField&>
+            fld.boundaryField();
-            (
-                fld.boundaryField()
-            );
 
         if (patchFieldDict.found(fld.name()))
         {
@@ -137,20 +131,17 @@ void Foam::fvMeshTools::setPatchFields
     const typename GeoField::value_type& value
 )
 {
-    HashTable<const GeoField*> flds
+    HashTable<GeoField*> flds
     (
         mesh.objectRegistry::lookupClass<GeoField>()
     );
 
-    forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
+    forAllIter(typename HashTable<GeoField*>, flds, iter)
     {
-        const GeoField& fld = *iter();
+        GeoField& fld = *iter();
 
         typename GeoField::GeometricBoundaryField& bfld =
-            const_cast<typename GeoField::GeometricBoundaryField&>
+            fld.boundaryField();
-            (
-                fld.boundaryField()
-            );
 
         bfld[patchI] == value;
     }
@@ -161,19 +152,15 @@ void Foam::fvMeshTools::setPatchFields
 template<class GeoField>
 void Foam::fvMeshTools::trimPatchFields(fvMesh& mesh, const label nPatches)
 {
-    HashTable<const GeoField*> flds
+    HashTable<GeoField*> flds
     (
         mesh.objectRegistry::lookupClass<GeoField>()
     );
 
-    forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
+    forAllIter(typename HashTable<GeoField*>, flds, iter)
     {
-        const GeoField& fld = *iter();
+        GeoField& fld = *iter();
-
+        fld.boundaryField().setSize(nPatches);
-        const_cast<typename GeoField::GeometricBoundaryField&>
-        (
-            fld.boundaryField()
-        ).setSize(nPatches);
     }
 }
 
@@ -186,20 +173,18 @@ void Foam::fvMeshTools::reorderPatchFields
     const labelList& oldToNew
 )
 {
-    HashTable<const GeoField*> flds
+    HashTable<GeoField*> flds
     (
         mesh.objectRegistry::lookupClass<GeoField>()
     );
 
-    forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
+    forAllIter(typename HashTable<GeoField*>, flds, iter)
     {
-        const GeoField& fld = *iter();
+        GeoField& fld = *iter();
 
         typename GeoField::GeometricBoundaryField& bfld =
-            const_cast<typename GeoField::GeometricBoundaryField&>
+            fld.boundaryField();
-            (
+
-                fld.boundaryField()
-            );
         bfld.reorder(oldToNew);
     }
 }

src/finiteVolume/Make/files

@@ -340,8 +340,7 @@ $(gradSchemes)/gaussGrad/gaussGrads.C
 
 $(gradSchemes)/leastSquaresGrad/leastSquaresVectors.C
 $(gradSchemes)/leastSquaresGrad/leastSquaresGrads.C
-$(gradSchemes)/extendedLeastSquaresGrad/extendedLeastSquaresVectors.C
+$(gradSchemes)/LeastSquaresGrad/LeastSquaresGrads.C
-$(gradSchemes)/extendedLeastSquaresGrad/extendedLeastSquaresGrads.C
 $(gradSchemes)/fourthGrad/fourthGrads.C
 
 limitedGradSchemes = $(gradSchemes)/limitedGradSchemes

src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControlTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,12 +34,12 @@ void Foam::solutionControl::storePrevIter() const
 {
     typedef GeometricField<Type, fvPatchField, volMesh> GeoField;
 
-    HashTable<const GeoField*>
+    HashTable<GeoField*>
         flds(mesh_.objectRegistry::lookupClass<GeoField>());
 
-    forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
+    forAllIter(typename HashTable<GeoField*>, flds, iter)
     {
-        GeoField& fld = const_cast<GeoField&>(*iter());
+        GeoField& fld = *iter();
 
         const word& fName = fld.name();
 

src/finiteVolume/finiteVolume/gradSchemes/LeastSquaresGrad/LeastSquaresGrad.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,18 +23,16 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "extendedLeastSquaresGrad.H"
+#include "LeastSquaresGrad.H"
-#include "extendedLeastSquaresVectors.H"
+#include "LeastSquaresVectors.H"
 #include "gaussGrad.H"
 #include "fvMesh.H"
 #include "volMesh.H"
-#include "surfaceMesh.H"
-#include "GeometricField.H"
 #include "zeroGradientFvPatchField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-template<class Type>
+template<class Type, class Stencil>
 Foam::tmp
 <
     Foam::GeometricField
@@ -44,15 +42,21 @@ Foam::tmp
     Foam::volMesh
     >
 >
-Foam::fv::extendedLeastSquaresGrad<Type>::calcGrad
+Foam::fv::LeastSquaresGrad<Type, Stencil>::calcGrad
 (
-    const GeometricField<Type, fvPatchField, volMesh>& vsf,
+    const GeometricField<Type, fvPatchField, volMesh>& vtf,
     const word& name
 ) const
 {
     typedef typename outerProduct<vector, Type>::type GradType;
 
-    const fvMesh& mesh = vsf.mesh();
+    const fvMesh& mesh = vtf.mesh();
+
+    // Get reference to least square vectors
+    const LeastSquaresVectors<Stencil>& lsv = LeastSquaresVectors<Stencil>::New
+    (
+        mesh
+    );
 
     tmp<GeometricField<GradType, fvPatchField, volMesh> > tlsGrad
     (
@@ -61,7 +65,7 @@ Foam::fv::extendedLeastSquaresGrad<Type>::calcGrad
             IOobject
             (
                 name,
-                vsf.instance(),
+                vtf.instance(),
                 mesh,
                 IOobject::NO_READ,
                 IOobject::NO_WRITE
@@ -70,87 +74,64 @@ Foam::fv::extendedLeastSquaresGrad<Type>::calcGrad
             dimensioned<GradType>
             (
                 "zero",
-                vsf.dimensions()/dimLength,
+                vtf.dimensions()/dimLength,
                 pTraits<GradType>::zero
             ),
             zeroGradientFvPatchField<GradType>::typeName
         )
     );
     GeometricField<GradType, fvPatchField, volMesh>& lsGrad = tlsGrad();
+    Field<GradType>& lsGradIf = lsGrad;
 
-    // Get reference to least square vectors
+    const extendedCentredCellToCellStencil& stencil = lsv.stencil();
-    const extendedLeastSquaresVectors& lsv = extendedLeastSquaresVectors::New
+    const List<List<label> >& stencilAddr = stencil.stencil();
-    (
+    const List<List<vector> >& lsvs = lsv.vectors();
-        mesh,
-        minDet_
-    );
 
-    const surfaceVectorField& ownLs = lsv.pVectors();
+    // Construct flat version of vtf
-    const surfaceVectorField& neiLs = lsv.nVectors();
+    // including all values referred to by the stencil
+    List<Type> flatVtf(stencil.map().constructSize(), pTraits<Type>::zero);
 
-    const labelUList& owner = mesh.owner();
+    // Insert internal values
-    const labelUList& neighbour = mesh.neighbour();
+    forAll(vtf, celli)
-
-    forAll(owner, facei)
     {
-        label own = owner[facei];
+        flatVtf[celli] = vtf[celli];
-        label nei = neighbour[facei];
-
-        Type deltaVsf = vsf[nei] - vsf[own];
-
-        lsGrad[own] += ownLs[facei]*deltaVsf;
-        lsGrad[nei] -= neiLs[facei]*deltaVsf;
     }
 
-    // Boundary faces
+    // Insert boundary values
-    forAll(vsf.boundaryField(), patchi)
+    forAll(vtf.boundaryField(), patchi)
     {
-        const fvsPatchVectorField& patchOwnLs = ownLs.boundaryField()[patchi];
+        const fvPatchField<Type>& ptf = vtf.boundaryField()[patchi];
 
-        const labelUList& faceCells =
+        label nCompact =
-            lsGrad.boundaryField()[patchi].patch().faceCells();
+            ptf.patch().start()
+          - mesh.nInternalFaces()
+          + mesh.nCells();
 
-        if (vsf.boundaryField()[patchi].coupled())
+        forAll(ptf, i)
         {
-            const Field<Type> neiVsf
+            flatVtf[nCompact++] = ptf[i];
-            (
-                vsf.boundaryField()[patchi].patchNeighbourField()
-            );
-
-            forAll(neiVsf, patchFaceI)
-            {
-                lsGrad[faceCells[patchFaceI]] +=
-                    patchOwnLs[patchFaceI]
-                   *(neiVsf[patchFaceI] - vsf[faceCells[patchFaceI]]);
         }
     }
-        else
+
+    // Do all swapping to complete flatVtf
+    stencil.map().distribute(flatVtf);
+
+    // Accumulate the cell-centred gradient from the
+    // weighted least-squares vectors and the flattened field values
+    forAll(stencilAddr, celli)
     {
-            const fvPatchField<Type>& patchVsf = vsf.boundaryField()[patchi];
+        const labelList& compactCells = stencilAddr[celli];
+        const List<vector>& lsvc = lsvs[celli];
 
-            forAll(patchVsf, patchFaceI)
+        forAll(compactCells, i)
         {
-                lsGrad[faceCells[patchFaceI]] +=
+            lsGradIf[celli] += lsvc[i]*flatVtf[compactCells[i]];
-                     patchOwnLs[patchFaceI]
-                    *(patchVsf[patchFaceI] - vsf[faceCells[patchFaceI]]);
-            }
         }
     }
 
-
+    // Correct the boundary conditions
-    const List<labelPair>& additionalCells = lsv.additionalCells();
-    const vectorField& additionalVectors = lsv.additionalVectors();
-
-    forAll(additionalCells, i)
-    {
-        lsGrad[additionalCells[i][0]] +=
-            additionalVectors[i]
-           *(vsf[additionalCells[i][1]] - vsf[additionalCells[i][0]]);
-    }
-
-
     lsGrad.correctBoundaryConditions();
-    gaussGrad<Type>::correctBoundaryConditions(vsf, lsGrad);
+    gaussGrad<Type>::correctBoundaryConditions(vtf, lsGrad);
 
     return tlsGrad;
 }

src/finiteVolume/finiteVolume/gradSchemes/LeastSquaresGrad/LeastSquaresGrad.H

@@ -0,0 +1,174 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::fv::LeastSquaresGrad
+
+Description
+    Gradient calculated using weighted least-squares on an arbitrary stencil.
+    The stencil type is provided via a template argument and any cell-based
+    stencil is supported:
+
+    \table
+        Stencil                     | Connections     | Scheme name
+        centredCFCCellToCellStencil | cell-face-cell  | Not Instantiated
+        centredCPCCellToCellStencil | cell-point-cell | pointCellsLeastSquares
+        centredCECCellToCellStencil | cell-edge-cell  | edgeCellsLeastSquares
+    \endtable
+
+    The first of these is not instantiated by default as the standard
+    leastSquaresGrad is equivalent and more efficient.
+
+    \heading Usage
+
+    Example of the gradient specification:
+    \verbatim
+    gradSchemes
+    {
+        default         pointCellsLeastSquares;
+    }
+    \endverbatim
+
+See Also
+    Foam::fv::LeastSquaresVectors
+    Foam::fv::leastSquaresGrad
+
+SourceFiles
+    LeastSquaresGrad.C
+    LeastSquaresVectors.H
+    LeastSquaresVectors.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef LeastSquaresGrad_H
+#define LeastSquaresGrad_H
+
+#include "gradScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class LeastSquaresGrad Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type, class Stencil>
+class LeastSquaresGrad
+:
+    public fv::gradScheme<Type>
+{
+    // Private Data
+
+        //- Minimum determinant criterion to choose extra cells
+        scalar minDet_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        LeastSquaresGrad(const LeastSquaresGrad&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const LeastSquaresGrad&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("LeastSquares");
+
+
+    // Constructors
+
+        //- Construct from Istream
+        LeastSquaresGrad(const fvMesh& mesh, Istream& schemeData)
+        :
+            gradScheme<Type>(mesh)
+        {}
+
+
+    // Member Functions
+
+        //- Return the gradient of the given field to the gradScheme::grad
+        //  for optional caching
+        virtual tmp
+        <
+            GeometricField
+            <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
+        > calcGrad
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& vsf,
+            const word& name
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Add the patch constructor functions to the hash tables
+
+#define makeLeastSquaresGradTypeScheme(SS, STENCIL, TYPE)                     \
+                                                                              \
+typedef LeastSquaresGrad<TYPE, STENCIL> LeastSquaresGrad##TYPE##STENCIL##_;   \
+defineTemplateTypeNameAndDebugWithName                                        \
+    (LeastSquaresGrad##TYPE##STENCIL##_, #SS, 0);                             \
+                                                                              \
+gradScheme<TYPE>::addIstreamConstructorToTable                                \
+<LeastSquaresGrad<TYPE, STENCIL> >                                            \
+    add##SS##STENCIL##TYPE##IstreamConstructorToTable_;
+
+#define makeLeastSquaresGradScheme(SS, STENCIL)                               \
+                                                                              \
+typedef LeastSquaresVectors<STENCIL> LeastSquaresVectors##STENCIL##_;         \
+defineTemplateTypeNameAndDebugWithName                                        \
+    (LeastSquaresVectors##STENCIL##_, #SS, 0);                                \
+                                                                              \
+makeLeastSquaresGradTypeScheme(SS,STENCIL,scalar)                             \
+makeLeastSquaresGradTypeScheme(SS,STENCIL,vector)
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "LeastSquaresGrad.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/finiteVolume/gradSchemes/LeastSquaresGrad/LeastSquaresGrads.C

@@ -0,0 +1,57 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "LeastSquaresGrad.H"
+//#include "centredCFCCellToCellStencilObject.H"
+#include "centredCPCCellToCellStencilObject.H"
+#include "centredCECCellToCellStencilObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+    // makeLeastSquaresGradScheme
+    // (
+    //     faceCellsLeastSquares,
+    //     centredCFCCellToCellStencilObject
+    // )
+
+    makeLeastSquaresGradScheme
+    (
+        pointCellsLeastSquares,
+        centredCPCCellToCellStencilObject
+    )
+
+    makeLeastSquaresGradScheme
+    (
+        edgeCellsLeastSquares,
+        centredCECCellToCellStencilObject
+    )
+}
+}
+
+// ************************************************************************* //

src/finiteVolume/finiteVolume/gradSchemes/LeastSquaresGrad/LeastSquaresVectors.C

@@ -0,0 +1,118 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "LeastSquaresVectors.H"
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+template<class Stencil>
+Foam::fv::LeastSquaresVectors<Stencil>::LeastSquaresVectors
+(
+    const fvMesh& mesh
+)
+:
+    MeshObject<fvMesh, Foam::MoveableMeshObject, LeastSquaresVectors>(mesh),
+    vectors_(mesh.nCells())
+{
+    calcLeastSquaresVectors();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+template<class Stencil>
+Foam::fv::LeastSquaresVectors<Stencil>::~LeastSquaresVectors()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Stencil>
+void Foam::fv::LeastSquaresVectors<Stencil>::calcLeastSquaresVectors()
+{
+    if (debug)
+    {
+        Info<< "LeastSquaresVectors::calcLeastSquaresVectors() :"
+            << "Calculating least square gradient vectors"
+            << endl;
+    }
+
+    const fvMesh& mesh = this->mesh_;
+    const extendedCentredCellToCellStencil& stencil = this->stencil();
+
+    stencil.collectData(mesh.C(), vectors_);
+
+    // Create the base form of the dd-tensor
+    // including components for the "empty" directions
+    symmTensor dd0(sqr((Vector<label>::one - mesh.geometricD())/2));
+
+    forAll (vectors_, i)
+    {
+        List<vector>& lsvi = vectors_[i];
+        symmTensor dd(dd0);
+
+        // The current cell is 0 in the stencil
+        // Calculate the deltas and sum the weighted dd
+        for (label j=1; j<lsvi.size(); j++)
+        {
+            lsvi[j] = lsvi[j] - lsvi[0];
+            scalar magSqrLsvi = magSqr(lsvi[j]);
+            dd += sqr(lsvi[j])/magSqrLsvi;
+            lsvi[j] /= magSqrLsvi;
+        }
+
+        // Invert dd
+        dd = inv(dd);
+
+        // Remove the components corresponding to the empty directions
+        dd -= dd0;
+
+        // Finalize the gradient weighting vectors
+        lsvi[0] = vector::zero;
+        for (label j=1; j<lsvi.size(); j++)
+        {
+            lsvi[j] = dd & lsvi[j];
+            lsvi[0] -= lsvi[j];
+        }
+    }
+
+    if (debug)
+    {
+        Info<< "LeastSquaresVectors::calcLeastSquaresVectors() :"
+            << "Finished calculating least square gradient vectors"
+            << endl;
+    }
+}
+
+
+template<class Stencil>
+bool Foam::fv::LeastSquaresVectors<Stencil>::movePoints()
+{
+    calcLeastSquaresVectors();
+    return true;
+}
+
+
+// ************************************************************************* //

src/finiteVolume/finiteVolume/gradSchemes/LeastSquaresGrad/LeastSquaresVectors.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,103 +22,110 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::extendedLeastSquaresVectors
+    Foam::fv::LeastSquaresVectors
 
 Description
-    Extended molecule least-squares gradient scheme vectors
+    Least-squares gradient scheme vectors
+
+See Also
+    Foam::fv::LeastSquaresGrad
 
 SourceFiles
-    extendedLeastSquaresVectors.C
+    LeastSquaresVectors.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef extendedLeastSquaresVectors_H
+#ifndef LeastSquaresVectors_H
-#define extendedLeastSquaresVectors_H
+#define LeastSquaresVectors_H
 
+#include "extendedCentredCellToCellStencil.H"
 #include "MeshObject.H"
-#include "fvMesh.H"
-#include "surfaceFieldsFwd.H"
-#include "labelPair.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fv
+{
+
 /*---------------------------------------------------------------------------*\
-                     Class extendedLeastSquaresVectors Declaration
+                     Class LeastSquaresVectors Declaration
 \*---------------------------------------------------------------------------*/
 
-class extendedLeastSquaresVectors
+template<class Stencil>
+class LeastSquaresVectors
 :
-    public MeshObject<fvMesh, extendedLeastSquaresVectors>
+    public MeshObject<fvMesh, MoveableMeshObject, LeastSquaresVectors<Stencil> >
 {
     // Private data
 
-        //- Minimum determinant criterion to choose extra cells
-        scalar minDet_;
-
         //- Least-squares gradient vectors
-        mutable surfaceVectorField* pVectorsPtr_;
+        List<List<vector> > vectors_;
-        mutable surfaceVectorField* nVectorsPtr_;
-
-        mutable List<labelPair>* additionalCellsPtr_;
-        mutable vectorField* additionalVectorsPtr_;
 
 
     // Private Member Functions
 
-        //- Construct Least-squares gradient vectors
+        //- Calculate Least-squares gradient vectors
-        void makeLeastSquaresVectors() const;
+        void calcLeastSquaresVectors();
 
 
 public:
 
     // Declare name of the class and its debug switch
-    TypeName("extendedLeastSquaresVectors");
+    TypeName("LeastSquaresVectors");
 
 
     // Constructors
 
         //- Construct given an fvMesh and the minimum determinant criterion
-        extendedLeastSquaresVectors
+        LeastSquaresVectors
         (
-            const fvMesh&,
+            const fvMesh&
-            const scalar minDet
         );
 
 
     //- Destructor
-    virtual ~extendedLeastSquaresVectors();
+    virtual ~LeastSquaresVectors();
 
 
     // Member functions
 
-        //- Return reference to owner least square vectors
+        //- Return reference to the stencil
-        const surfaceVectorField& pVectors() const;
+        const extendedCentredCellToCellStencil& stencil() const
+        {
+            return Stencil::New(this->mesh_);
+        }
 
-        //- Return reference to neighbour least square vectors
+        //- Return reference to the least square vectors
-        const surfaceVectorField& nVectors() const;
+        const List<List<vector> >& vectors() const
+        {
+            return vectors_;
+        }
 
-
+        //- Update the least square vectors when the mesh moves
-        //- Return reference to additional cells
-        const List<labelPair>& additionalCells() const;
-
-        //- Return reference to additional least square vectors
-        const vectorField& additionalVectors() const;
-
-
-        //- Delete the least square vectors when the mesh moves
         virtual bool movePoints();
 };
 
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef NoRepository
+#   include "LeastSquaresVectors.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //

src/finiteVolume/finiteVolume/gradSchemes/extendedLeastSquaresGrad/extendedLeastSquaresGrad.H

@@ -1,136 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
-     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::fv::extendedLeastSquaresGrad
-
-Description
-    Second-order gradient scheme using least-squares.
-
-SourceFiles
-    extendedLeastSquaresGrad.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef extendedLeastSquaresGrad_H
-#define extendedLeastSquaresGrad_H
-
-#include "gradScheme.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace fv
-{
-
-/*---------------------------------------------------------------------------*\
-                      Class extendedLeastSquaresGrad Declaration
-\*---------------------------------------------------------------------------*/
-
-template<class Type>
-class extendedLeastSquaresGrad
-:
-    public fv::gradScheme<Type>
-{
-    // Private Data
-
-        //- Minimum determinant criterion to choose extra cells
-        scalar minDet_;
-
-
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        extendedLeastSquaresGrad(const extendedLeastSquaresGrad&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const extendedLeastSquaresGrad&);
-
-
-public:
-
-    //- Runtime type information
-    TypeName("extendedLeastSquares");
-
-
-    // Constructors
-
-        //- Construct from Istream
-        extendedLeastSquaresGrad(const fvMesh& mesh, Istream& schemeData)
-        :
-            gradScheme<Type>(mesh),
-            minDet_(readScalar(schemeData))
-        {
-            if (minDet_ < 0)    //-for facearea weighted: || minDet_ > 8)
-            {
-                FatalIOErrorIn
-                (
-                    "extendedLeastSquaresGrad"
-                    "(const fvMesh&, Istream& schemeData)",
-                    schemeData
-                )   << "Minimum determinant = " << minDet_
-                    << " should be >= 0"    // and <= 8"
-                    << exit(FatalIOError);
-            }
-        }
-
-
-    // Member Functions
-
-        //- Return the gradient of the given field to the gradScheme::grad
-        //  for optional caching
-        virtual tmp
-        <
-            GeometricField
-            <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
-        > calcGrad
-        (
-            const GeometricField<Type, fvPatchField, volMesh>& vsf,
-            const word& name
-        ) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace fv
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-#   include "extendedLeastSquaresGrad.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

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

@@ -1,453 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
-     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "extendedLeastSquaresVectors.H"
-#include "surfaceFields.H"
-#include "volFields.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-defineTypeNameAndDebug(extendedLeastSquaresVectors, 0);
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
-
-Foam::extendedLeastSquaresVectors::extendedLeastSquaresVectors
-(
-    const fvMesh& mesh,
-    const scalar minDet
-)
-:
-    MeshObject<fvMesh, extendedLeastSquaresVectors>(mesh),
-    minDet_(minDet),
-    pVectorsPtr_(NULL),
-    nVectorsPtr_(NULL),
-    additionalCellsPtr_(NULL),
-    additionalVectorsPtr_(NULL)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
-
-Foam::extendedLeastSquaresVectors::~extendedLeastSquaresVectors()
-{
-    deleteDemandDrivenData(pVectorsPtr_);
-    deleteDemandDrivenData(nVectorsPtr_);
-
-    deleteDemandDrivenData(additionalCellsPtr_);
-    deleteDemandDrivenData(additionalVectorsPtr_);
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::extendedLeastSquaresVectors::makeLeastSquaresVectors() const
-{
-    if (debug)
-    {
-        Info<< "extendedLeastSquaresVectors::makeLeastSquaresVectors() :"
-            << "Constructing least square gradient vectors"
-            << endl;
-    }
-
-    const fvMesh& mesh = mesh_;
-
-    pVectorsPtr_ = new surfaceVectorField
-    (
-        IOobject
-        (
-            "LeastSquaresP",
-            mesh_.pointsInstance(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        mesh_,
-        dimensionedVector("zero", dimless/dimLength, vector::zero)
-    );
-    surfaceVectorField& lsP = *pVectorsPtr_;
-
-    nVectorsPtr_ = new surfaceVectorField
-    (
-        IOobject
-        (
-            "LeastSquaresN",
-            mesh_.pointsInstance(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        mesh_,
-        dimensionedVector("zero", dimless/dimLength, vector::zero)
-    );
-    surfaceVectorField& lsN = *nVectorsPtr_;
-
-    // Set local references to mesh data
-    const labelUList& owner = mesh_.owner();
-    const labelUList& neighbour = mesh_.neighbour();
-
-
-    // Determine number of dimensions and (for 2D) missing dimension
-    label nDims = 0;
-    label twoD = -1;
-    for (direction dir = 0; dir < vector::nComponents; dir++)
-    {
-        if (mesh.geometricD()[dir] == 1)
-        {
-            nDims++;
-        }
-        else
-        {
-            twoD = dir;
-        }
-    }
-
-    if (nDims == 1)
-    {
-        FatalErrorIn
-        (
-            "extendedLeastSquaresVectors::makeLeastSquaresVectors() const"
-        )   << "Found a mesh with only one geometric dimension : "
-            << mesh.geometricD()
-            << exit(FatalError);
-    }
-    else if (nDims == 2)
-    {
-        Info<< "extendedLeastSquares : detected " << nDims
-            << " valid directions. Missing direction " << twoD << endl;
-    }
-
-
-    const volVectorField& C = mesh.C();
-
-    // Set up temporary storage for the dd tensor (before inversion)
-    symmTensorField dd(mesh_.nCells(), symmTensor::zero);
-
-    forAll(owner, facei)
-    {
-        label own = owner[facei];
-        label nei = neighbour[facei];
-
-        vector d = C[nei] - C[own];
-
-        const symmTensor wdd(1.0/magSqr(d)*sqr(d));
-
-        dd[own] += wdd;
-        dd[nei] += wdd;
-    }
-
-    // Visit the boundaries. Coupled boundaries are taken into account
-    // in the construction of d vectors.
-    surfaceVectorField::GeometricBoundaryField& blsP = lsP.boundaryField();
-
-    forAll(blsP, patchi)
-    {
-        const fvPatch& p = blsP[patchi].patch();
-        const labelUList& faceCells = p.faceCells();
-
-        // Build the d-vectors
-        vectorField pd(p.delta());
-
-        forAll(pd, patchFaceI)
-        {
-            dd[faceCells[patchFaceI]] +=
-                (1.0/magSqr(pd[patchFaceI]))*sqr(pd[patchFaceI]);
-        }
-    }
-
-
-    // Check for missing dimensions
-    // Add the missing eigenvector (such that it does not
-    // affect the determinant)
-    if (nDims == 2)
-    {
-        forAll(dd, cellI)
-        {
-            if (twoD == 0)
-            {
-                dd[cellI].xx() = 1;
-            }
-            else if (twoD == 1)
-            {
-                dd[cellI].yy() = 1;
-            }
-            else
-            {
-                dd[cellI].zz() = 1;
-            }
-        }
-    }
-    scalarField detdd(det(dd));
-
-    Info<< "max(detdd) = " << max(detdd) << nl
-        << "min(detdd) = " << min(detdd) << nl
-        << "average(detdd) = " << average(detdd) << endl;
-
-    label nAdaptedCells = 0;
-    label nAddCells = 0;
-    label maxNaddCells = 4*detdd.size();
-    additionalCellsPtr_ = new List<labelPair>(maxNaddCells);
-    List<labelPair>& additionalCells_ = *additionalCellsPtr_;
-
-    forAll(detdd, i)
-    {
-        label count = 0;
-
-        label oldNAddCells = nAddCells;
-
-        while (++count < 100 && detdd[i] < minDet_)
-        {
-            if (nAddCells == maxNaddCells)
-            {
-                FatalErrorIn
-                (
-                    "extendedLeastSquaresVectors::"
-                    "makeLeastSquaresVectors() const"
-                )   << "nAddCells exceeds maxNaddCells ("
-                    << maxNaddCells << ")"
-                    << exit(FatalError);
-            }
-
-            labelList pointLabels = mesh.cells()[i].labels(mesh.faces());
-
-            scalar maxDetddij = 0.0;
-
-            label addCell = -1;
-
-            forAll(pointLabels, j)
-            {
-                forAll(mesh.pointCells()[pointLabels[j]], k)
-                {
-                    label cellj = mesh.pointCells()[pointLabels[j]][k];
-
-                    if (cellj != i)
-                    {
-                        vector dCij = (mesh.C()[cellj] - mesh.C()[i]);
-
-                        symmTensor ddij =
-                            dd[i] + (1.0/magSqr(dCij))*sqr(dCij);
-
-                        // Add the missing eigenvector (such that it does not
-                        // affect the determinant)
-                        if (nDims == 2)
-                        {
-                            if (twoD == 0)
-                            {
-                                ddij.xx() = 1;
-                            }
-                            else if (twoD == 1)
-                            {
-                                ddij.yy() = 1;
-                            }
-                            else
-                            {
-                                ddij.zz() = 1;
-                            }
-                        }
-
-                        scalar detddij = det(ddij);
-
-                        if (detddij > maxDetddij)
-                        {
-                            addCell = cellj;
-                            maxDetddij = detddij;
-                        }
-                    }
-                }
-            }
-
-            if (addCell != -1)
-            {
-                additionalCells_[nAddCells][0] = i;
-                additionalCells_[nAddCells++][1] = addCell;
-                vector dCij = mesh.C()[addCell] - mesh.C()[i];
-                dd[i] += (1.0/magSqr(dCij))*sqr(dCij);
-
-                // Add the missing eigenvector (such that it does not
-                // affect the determinant)
-                if (nDims == 2)
-                {
-                    if (twoD == 0)
-                    {
-                        dd[i].xx() = 1;
-                    }
-                    else if (twoD == 1)
-                    {
-                        dd[i].yy() = 1;
-                    }
-                    else
-                    {
-                        dd[i].zz() = 1;
-                    }
-                }
-
-                detdd[i] = det(dd[i]);
-            }
-        }
-
-        if (oldNAddCells < nAddCells)
-        {
-            nAdaptedCells++;
-        }
-    }
-
-    additionalCells_.setSize(nAddCells);
-
-    reduce(nAddCells, sumOp<label>());
-    reduce(nAdaptedCells, sumOp<label>());
-    if (nAddCells)
-    {
-        Info<< "max(detdd) = " << max(detdd) << nl
-            << "min(detdd) = " << min(detdd) << nl
-            << "average(detdd) = " << average(detdd) << nl
-            << "nAdapted/nCells = "
-            << scalar(nAdaptedCells)/mesh.globalData().nTotalCells() << nl
-            << "nAddCells/nCells = "
-            << scalar(nAddCells)/mesh.globalData().nTotalCells()
-            << endl;
-    }
-
-    // Invert the dd tensor
-    const symmTensorField invDd(inv(dd));
-
-
-    // Revisit all faces and calculate the lsP and lsN vectors
-    forAll(owner, facei)
-    {
-        label own = owner[facei];
-        label nei = neighbour[facei];
-
-        vector d = C[nei] - C[own];
-
-        lsP[facei] = (1.0/magSqr(d))*(invDd[own] & d);
-        lsN[facei] = ((-1.0)/magSqr(d))*(invDd[nei] & d);
-    }
-
-    forAll(blsP, patchI)
-    {
-        fvsPatchVectorField& patchLsP = blsP[patchI];
-
-        const fvPatch& p = patchLsP.patch();
-        const labelUList& faceCells = p.faceCells();
-
-        // Build the d-vectors
-        vectorField pd(p.delta());
-
-        forAll(p, patchFaceI)
-        {
-            patchLsP[patchFaceI] =
-                (1.0/magSqr(pd[patchFaceI]))
-               *(invDd[faceCells[patchFaceI]] & pd[patchFaceI]);
-        }
-    }
-
-
-    additionalVectorsPtr_ = new vectorField(additionalCells_.size());
-    vectorField& additionalVectors_ = *additionalVectorsPtr_;
-
-    forAll(additionalCells_, i)
-    {
-        vector dCij =
-            mesh.C()[additionalCells_[i][1]] - mesh.C()[additionalCells_[i][0]];
-
-        additionalVectors_[i] =
-            (1.0/magSqr(dCij))*(invDd[additionalCells_[i][0]] & dCij);
-    }
-
-    if (debug)
-    {
-        Info<< "extendedLeastSquaresVectors::makeLeastSquaresVectors() :"
-            << "Finished constructing least square gradient vectors"
-            << endl;
-    }
-}
-
-
-const Foam::surfaceVectorField&
-Foam::extendedLeastSquaresVectors::pVectors() const
-{
-    if (!pVectorsPtr_)
-    {
-        makeLeastSquaresVectors();
-    }
-
-    return *pVectorsPtr_;
-}
-
-
-const Foam::surfaceVectorField&
-Foam::extendedLeastSquaresVectors::nVectors() const
-{
-    if (!nVectorsPtr_)
-    {
-        makeLeastSquaresVectors();
-    }
-
-    return *nVectorsPtr_;
-}
-
-
-const Foam::List<Foam::labelPair>&
-Foam::extendedLeastSquaresVectors::additionalCells() const
-{
-    if (!additionalCellsPtr_)
-    {
-        makeLeastSquaresVectors();
-    }
-
-    return *additionalCellsPtr_;
-}
-
-
-const Foam::vectorField&
-Foam::extendedLeastSquaresVectors::additionalVectors() const
-{
-    if (!additionalVectorsPtr_)
-    {
-        makeLeastSquaresVectors();
-    }
-
-    return *additionalVectorsPtr_;
-}
-
-
-bool Foam::extendedLeastSquaresVectors::movePoints()
-{
-    deleteDemandDrivenData(pVectorsPtr_);
-    deleteDemandDrivenData(nVectorsPtr_);
-
-    deleteDemandDrivenData(additionalCellsPtr_);
-    deleteDemandDrivenData(additionalVectorsPtr_);
-
-    return true;
-}
-
-
-// ************************************************************************* //

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,14 +24,13 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "leastSquaresVectors.H"
-#include "surfaceFields.H"
 #include "volFields.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-defineTypeNameAndDebug(leastSquaresVectors, 0);
+    defineTypeNameAndDebug(leastSquaresVectors, 0);
 }
 
 
@@ -39,35 +38,8 @@ defineTypeNameAndDebug(leastSquaresVectors, 0);
 
 Foam::leastSquaresVectors::leastSquaresVectors(const fvMesh& mesh)
 :
-    MeshObject<fvMesh, leastSquaresVectors>(mesh),
+    MeshObject<fvMesh, Foam::MoveableMeshObject, leastSquaresVectors>(mesh),
-    pVectorsPtr_(NULL),
+    pVectors_
-    nVectorsPtr_(NULL)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
-
-Foam::leastSquaresVectors::~leastSquaresVectors()
-{
-    deleteDemandDrivenData(pVectorsPtr_);
-    deleteDemandDrivenData(nVectorsPtr_);
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
-{
-    if (debug)
-    {
-        Info<< "leastSquaresVectors::makeLeastSquaresVectors() :"
-            << "Constructing least square gradient vectors"
-            << endl;
-    }
-
-    const fvMesh& mesh = mesh_;
-
-    pVectorsPtr_ = new surfaceVectorField
     (
         IOobject
         (
@@ -80,10 +52,8 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
         ),
         mesh_,
         dimensionedVector("zero", dimless/dimLength, vector::zero)
-    );
+    ),
-    surfaceVectorField& lsP = *pVectorsPtr_;
+    nVectors_
-
-    nVectorsPtr_ = new surfaceVectorField
     (
         IOobject
         (
@@ -96,8 +66,30 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
         ),
         mesh_,
         dimensionedVector("zero", dimless/dimLength, vector::zero)
-    );
+    )
-    surfaceVectorField& lsN = *nVectorsPtr_;
+{
+    calcLeastSquaresVectors();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::leastSquaresVectors::~leastSquaresVectors()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::leastSquaresVectors::calcLeastSquaresVectors()
+{
+    if (debug)
+    {
+        Info<< "leastSquaresVectors::calcLeastSquaresVectors() :"
+            << "Calculating least square gradient vectors"
+            << endl;
+    }
+
+    const fvMesh& mesh = mesh_;
 
     // Set local references to mesh data
     const labelUList& owner = mesh_.owner();
@@ -124,7 +116,8 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
     }
 
 
-    surfaceVectorField::GeometricBoundaryField& blsP = lsP.boundaryField();
+    surfaceVectorField::GeometricBoundaryField& blsP =
+        pVectors_.boundaryField();
 
     forAll(blsP, patchi)
     {
@@ -164,7 +157,7 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
     const symmTensorField invDd(inv(dd));
 
 
-    // Revisit all faces and calculate the lsP and lsN vectors
+    // Revisit all faces and calculate the pVectors_ and nVectors_ vectors
     forAll(owner, facei)
     {
         label own = owner[facei];
@@ -173,8 +166,8 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
         vector d = C[nei] - C[own];
         scalar magSfByMagSqrd = magSf[facei]/magSqr(d);
 
-        lsP[facei] = (1 - w[facei])*magSfByMagSqrd*(invDd[own] & d);
+        pVectors_[facei] = (1 - w[facei])*magSfByMagSqrd*(invDd[own] & d);
-        lsN[facei] = -w[facei]*magSfByMagSqrd*(invDd[nei] & d);
+        nVectors_[facei] = -w[facei]*magSfByMagSqrd*(invDd[nei] & d);
     }
 
     forAll(blsP, patchi)
@@ -214,127 +207,18 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
         }
     }
 
-
-    // For 3D meshes check the determinant of the dd tensor and switch to
-    // Gauss if it is less than 3
-    /* Currently the det(dd[celli]) criterion is incorrect: dd is weighted by Sf
-    if (mesh.nGeometricD() == 3)
-    {
-        label nBadCells = 0;
-
-        const cellList& cells = mesh.cells();
-        const scalarField& V = mesh.V();
-        const surfaceVectorField& Sf = mesh.Sf();
-        const surfaceScalarField& w = mesh.weights();
-
-        forAll(dd, celli)
-        {
-            if (det(dd[celli]) < 3)
-            {
-                nBadCells++;
-
-                const cell& c = cells[celli];
-
-                forAll(c, cellFacei)
-                {
-                    label facei = c[cellFacei];
-
-                    if (mesh.isInternalFace(facei))
-                    {
-                        scalar wf = max(min(w[facei], 0.8), 0.2);
-
-                        if (celli == owner[facei])
-                        {
-                            lsP[facei] = (1 - wf)*Sf[facei]/V[celli];
-                        }
-                        else
-                        {
-                            lsN[facei] = -wf*Sf[facei]/V[celli];
-                        }
-                    }
-                    else
-                    {
-                        label patchi = mesh.boundaryMesh().whichPatch(facei);
-
-                        if (mesh.boundary()[patchi].size())
-                        {
-                            label patchFacei =
-                                facei - mesh.boundaryMesh()[patchi].start();
-
-                            if (w.boundaryField()[patchi].coupled())
-                            {
-                                scalar wf = max
-                                (
-                                    min
-                                    (
-                                        w.boundaryField()[patchi][patchFacei],
-                                        0.8
-                                    ),
-                                    0.2
-                                );
-
-                                lsP.boundaryField()[patchi][patchFacei] =
-                                    (1 - wf)
-                                   *Sf.boundaryField()[patchi][patchFacei]
-                                   /V[celli];
-                            }
-                            else
-                            {
-                                lsP.boundaryField()[patchi][patchFacei] =
-                                    Sf.boundaryField()[patchi][patchFacei]
-                                   /V[celli];
-                            }
-                        }
-                    }
-                }
-            }
-        }
-
     if (debug)
     {
-            InfoIn("leastSquaresVectors::makeLeastSquaresVectors()")
+        Info<< "leastSquaresVectors::calcLeastSquaresVectors() :"
-                << "number of bad cells switched to Gauss = " << nBadCells
+            << "Finished calculating least square gradient vectors"
             << endl;
     }
-    }
-    */
-
-    if (debug)
-    {
-        Info<< "leastSquaresVectors::makeLeastSquaresVectors() :"
-            << "Finished constructing least square gradient vectors"
-            << endl;
-    }
-}
-
-
-const Foam::surfaceVectorField& Foam::leastSquaresVectors::pVectors() const
-{
-    if (!pVectorsPtr_)
-    {
-        makeLeastSquaresVectors();
-    }
-
-    return *pVectorsPtr_;
-}
-
-
-const Foam::surfaceVectorField& Foam::leastSquaresVectors::nVectors() const
-{
-    if (!nVectorsPtr_)
-    {
-        makeLeastSquaresVectors();
-    }
-
-    return *nVectorsPtr_;
 }
 
 
 bool Foam::leastSquaresVectors::movePoints()
 {
-    deleteDemandDrivenData(pVectorsPtr_);
+    calcLeastSquaresVectors();
-    deleteDemandDrivenData(nVectorsPtr_);
-
     return true;
 }
 

src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,8 +37,7 @@ SourceFiles
 
 #include "MeshObject.H"
 #include "fvMesh.H"
-#include "surfaceFieldsFwd.H"
+#include "surfaceFields.H"
-#include "labelPair.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,19 +50,19 @@ namespace Foam
 
 class leastSquaresVectors
 :
-    public MeshObject<fvMesh, leastSquaresVectors>
+    public MeshObject<fvMesh, MoveableMeshObject, leastSquaresVectors>
 {
     // Private data
 
         //- Least-squares gradient vectors
-        mutable surfaceVectorField* pVectorsPtr_;
+        surfaceVectorField pVectors_;
-        mutable surfaceVectorField* nVectorsPtr_;
+        surfaceVectorField nVectors_;
 
 
     // Private Member Functions
 
         //- Construct Least-squares gradient vectors
-        void makeLeastSquaresVectors() const;
+        void calcLeastSquaresVectors();
 
 
 public:
@@ -85,11 +84,16 @@ public:
     // Member functions
 
         //- Return reference to owner least square vectors
-        const surfaceVectorField& pVectors() const;
+        const surfaceVectorField& pVectors() const
+        {
+            return pVectors_;
+        }
 
         //- Return reference to neighbour least square vectors
-        const surfaceVectorField& nVectors() const;
+        const surfaceVectorField& nVectors() const
-
+        {
+            return nVectors_;
+        }
 
         //- Delete the least square vectors when the mesh moves
         virtual bool movePoints();

src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -601,12 +601,13 @@ void Foam::MULES::limiter
         {
             fvsPatchScalarField& lambdaPf = lambdaBf[patchi];
             const scalarField& phiCorrfPf = phiCorrBf[patchi];
+            const fvPatchScalarField& psiPf = psiBf[patchi];
 
             if (isA<wedgeFvPatch>(mesh.boundary()[patchi]))
             {
                 lambdaPf = 0;
             }
-            else
+            else if (psiPf.coupled())
             {
                 const labelList& pFaceCells =
                     mesh.boundary()[patchi].faceCells();
@@ -627,6 +628,32 @@ void Foam::MULES::limiter
                     }
                 }
             }
+            else
+            {
+                const labelList& pFaceCells =
+                    mesh.boundary()[patchi].faceCells();
+                const scalarField& phiBDPf = phiBDBf[patchi];
+
+                forAll(lambdaPf, pFacei)
+                {
+                    // Limit outlet faces only
+                    if (phiBDPf[pFacei] > 0)
+                    {
+                        label pfCelli = pFaceCells[pFacei];
+
+                        if (phiCorrfPf[pFacei] > 0.0)
+                        {
+                            lambdaPf[pFacei] =
+                                min(lambdaPf[pFacei], lambdap[pfCelli]);
+                        }
+                        else
+                        {
+                            lambdaPf[pFacei] =
+                                min(lambdaPf[pFacei], lambdam[pfCelli]);
+                        }
+                    }
+                }
+            }
         }
 
         syncTools::syncFaceList(mesh, allLambda, minEqOp<scalar>());

src/finiteVolume/fvMesh/extendedStencil/cellToCell/MeshObjects/centredCECCellToCellStencilObject.H

@@ -48,7 +48,12 @@ namespace Foam
 
 class centredCECCellToCellStencilObject
 :
-    public MeshObject<fvMesh, centredCECCellToCellStencilObject>,
+    public MeshObject
+    <
+        fvMesh,
+        TopologicalMeshObject,
+        centredCECCellToCellStencilObject
+    >,
     public extendedCentredCellToCellStencil
 {
 
@@ -64,7 +69,12 @@ public:
             const fvMesh& mesh
         )
         :
-            MeshObject<fvMesh, centredCECCellToCellStencilObject>(mesh),
+            MeshObject
+            <
+                fvMesh,
+                Foam::TopologicalMeshObject,
+                centredCECCellToCellStencilObject
+            >(mesh),
             extendedCentredCellToCellStencil(CECCellToCellStencil(mesh))
         {}
 

src/finiteVolume/fvMesh/extendedStencil/cellToCell/MeshObjects/centredCFCCellToCellStencilObject.H

@@ -48,7 +48,12 @@ namespace Foam
 
 class centredCFCCellToCellStencilObject
 :
-    public MeshObject<fvMesh, centredCFCCellToCellStencilObject>,
+    public MeshObject
+    <
+        fvMesh,
+        TopologicalMeshObject,
+        centredCFCCellToCellStencilObject
+    >,
     public extendedCentredCellToCellStencil
 {
 
@@ -64,7 +69,12 @@ public:
             const fvMesh& mesh
         )
         :
-            MeshObject<fvMesh, centredCFCCellToCellStencilObject>(mesh),
+            MeshObject
+            <
+                fvMesh,
+                Foam::TopologicalMeshObject,
+                centredCFCCellToCellStencilObject
+            >(mesh),
             extendedCentredCellToCellStencil(CFCCellToCellStencil(mesh))
         {}
 

src/finiteVolume/fvMesh/extendedStencil/cellToCell/MeshObjects/centredCPCCellToCellStencilObject.H

@@ -48,7 +48,12 @@ namespace Foam
 
 class centredCPCCellToCellStencilObject
 :
-    public MeshObject<fvMesh, centredCPCCellToCellStencilObject>,
+    public MeshObject
+    <
+        fvMesh,
+        TopologicalMeshObject,
+        centredCPCCellToCellStencilObject
+    >,
     public extendedCentredCellToCellStencil
 {
 
@@ -64,7 +69,12 @@ public:
             const fvMesh& mesh
         )
         :
-            MeshObject<fvMesh, centredCPCCellToCellStencilObject>(mesh),
+            MeshObject
+            <
+                fvMesh,
+                Foam::TopologicalMeshObject,
+                centredCPCCellToCellStencilObject
+            >(mesh),
             extendedCentredCellToCellStencil(CPCCellToCellStencil(mesh))
         {}
 

src/finiteVolume/fvMesh/extendedStencil/cellToCell/extendedCellToCellStencilTemplates.C

@@ -76,15 +76,14 @@ Foam::tmp
     );
     WeightedFieldType& wf = twf();
 
-    // cells
+    forAll(wf, celli)
-    forAll(wf, cellI)
     {
-        const List<Type>& stField = stencilFld[cellI];
+        const List<Type>& stField = stencilFld[celli];
-        const List<WeightType>& stWeight = stencilWeights[cellI];
+        const List<WeightType>& stWeight = stencilWeights[celli];
 
         forAll(stField, i)
         {
-            wf[cellI] += stWeight[i]*stField[i];
+            wf[celli] += stWeight[i]*stField[i];
         }
     }
 

src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/centredCECCellToFaceStencilObject.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -48,7 +48,12 @@ namespace Foam
 
 class centredCECCellToFaceStencilObject
 :
-    public MeshObject<fvMesh, centredCECCellToFaceStencilObject>,
+    public MeshObject
+    <
+        fvMesh,
+        TopologicalMeshObject,
+        centredCECCellToFaceStencilObject
+    >,
     public extendedCentredCellToFaceStencil
 {
 
@@ -64,7 +69,12 @@ public:
             const fvMesh& mesh
         )
         :
-            MeshObject<fvMesh, centredCECCellToFaceStencilObject>(mesh),
+            MeshObject
+            <
+                fvMesh,
+                Foam::TopologicalMeshObject,
+                centredCECCellToFaceStencilObject
+            >(mesh),
             extendedCentredCellToFaceStencil(CECCellToFaceStencil(mesh))
         {
             if (extendedCellToFaceStencil::debug)

src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/centredCFCCellToFaceStencilObject.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -48,7 +48,12 @@ namespace Foam
 
 class centredCFCCellToFaceStencilObject
 :
-    public MeshObject<fvMesh, centredCFCCellToFaceStencilObject>,
+    public MeshObject
+    <
+        fvMesh,
+        TopologicalMeshObject,
+        centredCFCCellToFaceStencilObject
+    >,
     public extendedCentredCellToFaceStencil
 {
 
@@ -64,7 +69,12 @@ public:
             const fvMesh& mesh
         )
         :
-            MeshObject<fvMesh, centredCFCCellToFaceStencilObject>(mesh),
+            MeshObject
+            <
+                fvMesh,
+                Foam::TopologicalMeshObject,
+                centredCFCCellToFaceStencilObject
+            >(mesh),
             extendedCentredCellToFaceStencil(CFCCellToFaceStencil(mesh))
         {
             if (extendedCellToFaceStencil::debug)

src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/centredCPCCellToFaceStencilObject.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -48,7 +48,12 @@ namespace Foam
 
 class centredCPCCellToFaceStencilObject
 :
-    public MeshObject<fvMesh, centredCPCCellToFaceStencilObject>,
+    public MeshObject
+    <
+        fvMesh,
+        TopologicalMeshObject,
+        centredCPCCellToFaceStencilObject
+    >,
     public extendedCentredCellToFaceStencil
 {
 
@@ -64,7 +69,12 @@ public:
             const fvMesh& mesh
         )
         :
-            MeshObject<fvMesh, centredCPCCellToFaceStencilObject>(mesh),
+            MeshObject
+            <
+                fvMesh,
+                Foam::TopologicalMeshObject,
+                centredCPCCellToFaceStencilObject
+            >(mesh),
             extendedCentredCellToFaceStencil(CPCCellToFaceStencil(mesh))
         {
             if (extendedCellToFaceStencil::debug)

src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/centredFECCellToFaceStencilObject.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -48,7 +48,12 @@ namespace Foam
 
 class centredFECCellToFaceStencilObject
 :
-    public MeshObject<fvMesh, centredFECCellToFaceStencilObject>,
+    public MeshObject
+    <
+        fvMesh,
+        TopologicalMeshObject,
+        centredFECCellToFaceStencilObject
+    >,
     public extendedCentredCellToFaceStencil
 {
 
@@ -64,7 +69,12 @@ public:
             const fvMesh& mesh
         )
         :
-            MeshObject<fvMesh, centredFECCellToFaceStencilObject>(mesh),
+            MeshObject
+            <
+                fvMesh,
+                Foam::TopologicalMeshObject,
+                centredFECCellToFaceStencilObject
+            >(mesh),
             extendedCentredCellToFaceStencil(FECCellToFaceStencil(mesh))
         {
             if (extendedCellToFaceStencil::debug)

src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/pureUpwindCFCCellToFaceStencilObject.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -48,7 +48,12 @@ namespace Foam
 
 class pureUpwindCFCCellToFaceStencilObject
 :
-    public MeshObject<fvMesh, pureUpwindCFCCellToFaceStencilObject>,
+    public MeshObject
+    <
+        fvMesh,
+        TopologicalMeshObject,
+        pureUpwindCFCCellToFaceStencilObject
+    >,
     public extendedUpwindCellToFaceStencil
 {
 
@@ -64,7 +69,12 @@ public:
             const fvMesh& mesh
         )
         :
-            MeshObject<fvMesh, pureUpwindCFCCellToFaceStencilObject>(mesh),
+            MeshObject
+            <
+                fvMesh,
+                Foam::TopologicalMeshObject,
+                pureUpwindCFCCellToFaceStencilObject
+            >(mesh),
             extendedUpwindCellToFaceStencil(CFCCellToFaceStencil(mesh))
         {
             if (extendedCellToFaceStencil::debug)

src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/upwindCECCellToFaceStencilObject.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -48,7 +48,12 @@ namespace Foam
 
 class upwindCECCellToFaceStencilObject
 :
-    public MeshObject<fvMesh, upwindCECCellToFaceStencilObject>,
+    public MeshObject
+    <
+        fvMesh,
+        TopologicalMeshObject,
+        upwindCECCellToFaceStencilObject
+    >,
     public extendedUpwindCellToFaceStencil
 {
 
@@ -66,7 +71,12 @@ public:
             const scalar minOpposedness
         )
         :
-            MeshObject<fvMesh, upwindCECCellToFaceStencilObject>(mesh),
+            MeshObject
+            <
+                fvMesh,
+                Foam::TopologicalMeshObject,
+                upwindCECCellToFaceStencilObject
+            >(mesh),
             extendedUpwindCellToFaceStencil
             (
                 CECCellToFaceStencil(mesh),

src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/upwindCFCCellToFaceStencilObject.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -48,7 +48,12 @@ namespace Foam
 
 class upwindCFCCellToFaceStencilObject
 :
-    public MeshObject<fvMesh, upwindCFCCellToFaceStencilObject>,
+    public MeshObject
+    <
+        fvMesh,
+        TopologicalMeshObject,
+        upwindCFCCellToFaceStencilObject
+    >,
     public extendedUpwindCellToFaceStencil
 {
 
@@ -66,7 +71,12 @@ public:
             const scalar minOpposedness
         )
         :
-            MeshObject<fvMesh, upwindCFCCellToFaceStencilObject>(mesh),
+            MeshObject
+            <
+                fvMesh,
+                Foam::TopologicalMeshObject,
+                upwindCFCCellToFaceStencilObject
+            >(mesh),
             extendedUpwindCellToFaceStencil
             (
                 CFCCellToFaceStencil(mesh),

src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/upwindCPCCellToFaceStencilObject.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -48,7 +48,12 @@ namespace Foam
 
 class upwindCPCCellToFaceStencilObject
 :
-    public MeshObject<fvMesh, upwindCPCCellToFaceStencilObject>,
+    public MeshObject
+    <
+        fvMesh,
+        TopologicalMeshObject,
+        upwindCPCCellToFaceStencilObject
+    >,
     public extendedUpwindCellToFaceStencil
 {
 
@@ -66,7 +71,12 @@ public:
             const scalar minOpposedness
         )
         :
-            MeshObject<fvMesh, upwindCPCCellToFaceStencilObject>(mesh),
+            MeshObject
+            <
+                fvMesh,
+                Foam::TopologicalMeshObject,
+                upwindCPCCellToFaceStencilObject
+            >(mesh),
             extendedUpwindCellToFaceStencil
             (
                 CPCCellToFaceStencil(mesh),

src/finiteVolume/fvMesh/extendedStencil/cellToFace/MeshObjects/upwindFECCellToFaceStencilObject.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -48,7 +48,12 @@ namespace Foam
 
 class upwindFECCellToFaceStencilObject
 :
-    public MeshObject<fvMesh, upwindFECCellToFaceStencilObject>,
+    public MeshObject
+    <
+        fvMesh,
+        TopologicalMeshObject,
+        upwindFECCellToFaceStencilObject
+    >,
     public extendedUpwindCellToFaceStencil
 {
 
@@ -66,7 +71,12 @@ public:
             const scalar minOpposedness
         )
         :
-            MeshObject<fvMesh, upwindFECCellToFaceStencilObject>(mesh),
+            MeshObject
+            <
+                fvMesh,
+                Foam::TopologicalMeshObject,
+                upwindFECCellToFaceStencilObject
+            >(mesh),
             extendedUpwindCellToFaceStencil
             (
                 FECCellToFaceStencil(mesh),

src/finiteVolume/fvMesh/extendedStencil/faceToCell/MeshObjects/centredCFCFaceToCellStencilObject.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -48,7 +48,12 @@ namespace Foam
 
 class centredCFCFaceToCellStencilObject
 :
-    public MeshObject<fvMesh, centredCFCFaceToCellStencilObject>,
+    public MeshObject
+    <
+        fvMesh,
+        TopologicalMeshObject,
+        centredCFCFaceToCellStencilObject
+    >,
     public extendedCentredFaceToCellStencil
 {
 
@@ -64,7 +69,12 @@ public:
             const fvMesh& mesh
         )
         :
-            MeshObject<fvMesh, centredCFCFaceToCellStencilObject>(mesh),
+            MeshObject
+            <
+                fvMesh,
+                Foam::TopologicalMeshObject,
+                centredCFCFaceToCellStencilObject
+            >(mesh),
             extendedCentredFaceToCellStencil(CFCFaceToCellStencil(mesh))
         {}
 

src/finiteVolume/fvMesh/fvMesh.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,42 +31,17 @@ License
 #include "SubField.H"
 #include "demandDrivenData.H"
 #include "fvMeshLduAddressing.H"
-#include "emptyPolyPatch.H"
 #include "mapPolyMesh.H"
 #include "MapFvFields.H"
 #include "fvMeshMapper.H"
 #include "mapClouds.H"
-
+#include "MeshObject.H"
-#include "volPointInterpolation.H"
-#include "extendedLeastSquaresVectors.H"
-#include "extendedLeastSquaresVectors.H"
-#include "leastSquaresVectors.H"
-#include "CentredFitData.H"
-#include "linearFitPolynomial.H"
-#include "quadraticFitPolynomial.H"
-#include "quadraticLinearFitPolynomial.H"
-//#include "quadraticFitSnGradData.H"
-#include "skewCorrectionVectors.H"
-
-
-#include "centredCECCellToFaceStencilObject.H"
-#include "centredCFCCellToFaceStencilObject.H"
-#include "centredCPCCellToFaceStencilObject.H"
-#include "centredFECCellToFaceStencilObject.H"
-#include "upwindCECCellToFaceStencilObject.H"
-#include "upwindCFCCellToFaceStencilObject.H"
-#include "upwindCPCCellToFaceStencilObject.H"
-#include "upwindFECCellToFaceStencilObject.H"
-
-#include "centredCFCFaceToCellStencilObject.H"
-#include "meshSearchMeshObject.H"
-#include "meshSearchFACECENTRETETSMeshObject.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-defineTypeNameAndDebug(fvMesh, 0);
+    defineTypeNameAndDebug(fvMesh, 0);
 }
 
 
@@ -74,6 +49,8 @@ defineTypeNameAndDebug(fvMesh, 0);
 
 void Foam::fvMesh::clearGeomNotOldVol()
 {
+    meshObject::clear<fvMesh, GeometricMeshObject>(*this);
+
     slicedVolScalarField::DimensionedInternalField* VPtr =
         static_cast<slicedVolScalarField::DimensionedInternalField*>(VPtr_);
     deleteDemandDrivenData(VPtr);
@@ -129,55 +106,13 @@ void Foam::fvMesh::clearGeom()
 
     // Mesh motion flux cannot be deleted here because the old-time flux
     // needs to be saved.
-
-    // Things geometry dependent that are not updated.
-    volPointInterpolation::Delete(*this);
-    extendedLeastSquaresVectors::Delete(*this);
-    leastSquaresVectors::Delete(*this);
-    CentredFitData<linearFitPolynomial>::Delete(*this);
-    CentredFitData<quadraticFitPolynomial>::Delete(*this);
-    CentredFitData<quadraticLinearFitPolynomial>::Delete(*this);
-    skewCorrectionVectors::Delete(*this);
-    //quadraticFitSnGradData::Delete(*this);
-
-    // Note: should be in polyMesh::clearGeom but meshSearch not in OpenFOAM
-    // library
-    meshSearchMeshObject::Delete(*this);
-    meshSearchFACECENTRETETSMeshObject::Delete(*this);
 }
 
 
 void Foam::fvMesh::clearAddressing()
 {
+    meshObject::clear<fvMesh, TopologicalMeshObject>(*this);
     deleteDemandDrivenData(lduPtr_);
-
-    // Hack until proper callbacks. Below are all the fvMesh-MeshObjects.
-
-    volPointInterpolation::Delete(*this);
-    extendedLeastSquaresVectors::Delete(*this);
-    leastSquaresVectors::Delete(*this);
-    CentredFitData<linearFitPolynomial>::Delete(*this);
-    CentredFitData<quadraticFitPolynomial>::Delete(*this);
-    CentredFitData<quadraticLinearFitPolynomial>::Delete(*this);
-    skewCorrectionVectors::Delete(*this);
-    //quadraticFitSnGradData::Delete(*this);
-
-    centredCECCellToFaceStencilObject::Delete(*this);
-    centredCFCCellToFaceStencilObject::Delete(*this);
-    centredCPCCellToFaceStencilObject::Delete(*this);
-    centredFECCellToFaceStencilObject::Delete(*this);
-    // Is this geometry related - cells distorting to upwind direction?
-    upwindCECCellToFaceStencilObject::Delete(*this);
-    upwindCFCCellToFaceStencilObject::Delete(*this);
-    upwindCPCCellToFaceStencilObject::Delete(*this);
-    upwindFECCellToFaceStencilObject::Delete(*this);
-
-    centredCFCFaceToCellStencilObject::Delete(*this);
-
-    // Note: should be in polyMesh::clearGeom but meshSearch not in OpenFOAM
-    // library
-    meshSearchMeshObject::Delete(*this);
-    meshSearchFACECENTRETETSMeshObject::Delete(*this);
 }
 
 
@@ -603,24 +538,6 @@ void Foam::fvMesh::mapFields(const mapPolyMesh& meshMap)
 }
 
 
-// Temporary helper function to call move points on
-// MeshObjects
-template<class Type>
-void MeshObjectMovePoints(const Foam::fvMesh& mesh)
-{
-    if (mesh.thisDb().foundObject<Type>(Type::typeName))
-    {
-        const_cast<Type&>
-        (
-            mesh.thisDb().lookupObject<Type>
-            (
-                Type::typeName
-            )
-        ).movePoints();
-    }
-}
-
-
 Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
 {
     // Grab old time volumes if the time has been incremented
@@ -714,17 +631,7 @@ Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
     boundary_.movePoints();
     surfaceInterpolation::movePoints();
 
-
+    meshObject::movePoints<fvMesh>(*this);
-    // Hack until proper callbacks. Below are all the fvMesh MeshObjects with a
-    // movePoints function.
-    MeshObjectMovePoints<volPointInterpolation>(*this);
-    MeshObjectMovePoints<extendedLeastSquaresVectors>(*this);
-    MeshObjectMovePoints<leastSquaresVectors>(*this);
-    MeshObjectMovePoints<CentredFitData<linearFitPolynomial> >(*this);
-    MeshObjectMovePoints<CentredFitData<quadraticFitPolynomial> >(*this);
-    MeshObjectMovePoints<CentredFitData<quadraticLinearFitPolynomial> >(*this);
-    MeshObjectMovePoints<skewCorrectionVectors>(*this);
-    //MeshObjectMovePoints<quadraticFitSnGradData>(*this);
 
     return tsweptVols;
 }
@@ -746,9 +653,7 @@ void Foam::fvMesh::updateMesh(const mapPolyMesh& mpm)
 
     clearAddressing();
 
-    // handleMorph() should also clear out the surfaceInterpolation.
+    meshObject::updateMesh<fvMesh>(*this, mpm);
-    // This is a temporary solution
-    surfaceInterpolation::movePoints();
 }
 
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -40,7 +40,7 @@ Foam::FitData<Form, ExtendedStencil, Polynomial>::FitData
     const scalar centralWeight
 )
 :
-    MeshObject<fvMesh, Form>(mesh),
+    MeshObject<fvMesh, Foam::MoveableMeshObject, Form>(mesh),
     stencil_(stencil),
     linearCorrection_(linearCorrection),
     linearLimitFactor_(linearLimitFactor),

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -54,7 +54,7 @@ namespace Foam
 template<class FitDataType, class ExtendedStencil, class Polynomial>
 class FitData
 :
-    public MeshObject<fvMesh, FitDataType>
+    public MeshObject<fvMesh, MoveableMeshObject, FitDataType>
 {
     // Private data
 

src/finiteVolume/interpolation/surfaceInterpolation/schemes/skewCorrected/skewCorrectionVectors.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,7 +24,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "skewCorrectionVectors.H"
-#include "surfaceFields.H"
 #include "volFields.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -39,28 +38,9 @@ namespace Foam
 
 Foam::skewCorrectionVectors::skewCorrectionVectors(const fvMesh& mesh)
 :
-    MeshObject<fvMesh, skewCorrectionVectors>(mesh),
+    MeshObject<fvMesh, Foam::MoveableMeshObject, skewCorrectionVectors>(mesh),
-    skew_(true),
+    skew_(false),
-    skewCorrectionVectors_(NULL)
+    skewCorrectionVectors_
-{}
-
-
-Foam::skewCorrectionVectors::~skewCorrectionVectors()
-{
-    deleteDemandDrivenData(skewCorrectionVectors_);
-}
-
-
-void Foam::skewCorrectionVectors::makeSkewCorrectionVectors() const
-{
-    if (debug)
-    {
-        Info<< "surfaceInterpolation::makeSkewCorrectionVectors() : "
-            << "Constructing skew correction vectors"
-            << endl;
-    }
-
-    skewCorrectionVectors_ = new surfaceVectorField
     (
         IOobject
         (
@@ -73,8 +53,24 @@ void Foam::skewCorrectionVectors::makeSkewCorrectionVectors() const
         ),
         mesh_,
         dimless
-    );
+    )
-    surfaceVectorField& SkewCorrVecs = *skewCorrectionVectors_;
+{
+    calcSkewCorrectionVectors();
+}
+
+
+Foam::skewCorrectionVectors::~skewCorrectionVectors()
+{}
+
+
+void Foam::skewCorrectionVectors::calcSkewCorrectionVectors()
+{
+    if (debug)
+    {
+        Info<< "surfaceInterpolation::calcSkewCorrectionVectors() : "
+            << "Calculating skew correction vectors"
+            << endl;
+    }
 
     // Set local references to mesh data
     const volVectorField& C = mesh_.C();
@@ -92,14 +88,15 @@ void Foam::skewCorrectionVectors::makeSkewCorrectionVectors() const
         vector d = C[nei] - C[own];
         vector Cpf = Cf[facei] - C[own];
 
-        SkewCorrVecs[facei] = Cpf - ((Sf[facei] & Cpf)/(Sf[facei] & d))*d;
+        skewCorrectionVectors_[facei] =
+            Cpf - ((Sf[facei] & Cpf)/(Sf[facei] & d))*d;
     }
 
 
-    forAll(SkewCorrVecs.boundaryField(), patchI)
+    forAll(skewCorrectionVectors_.boundaryField(), patchI)
     {
         fvsPatchVectorField& patchSkewCorrVecs =
-            SkewCorrVecs.boundaryField()[patchI];
+            skewCorrectionVectors_.boundaryField()[patchI];
 
         if (!patchSkewCorrVecs.coupled())
         {
@@ -132,21 +129,19 @@ void Foam::skewCorrectionVectors::makeSkewCorrectionVectors() const
 
     if (Sf.internalField().size())
     {
-        skewCoeff = max(mag(SkewCorrVecs)*mesh_.deltaCoeffs()).value();
+        skewCoeff =
+            max(mag(skewCorrectionVectors_)*mesh_.deltaCoeffs()).value();
     }
 
     if (debug)
     {
-        Info<< "surfaceInterpolation::makeSkewCorrectionVectors() : "
+        Info<< "surfaceInterpolation::calcSkewCorrectionVectors() : "
             << "skew coefficient = " << skewCoeff << endl;
     }
 
-    //skewCoeff = 0.0;
-
     if (skewCoeff < 1e-5)
     {
         skew_ = false;
-        deleteDemandDrivenData(skewCorrectionVectors_);
     }
     else
     {
@@ -155,43 +150,16 @@ void Foam::skewCorrectionVectors::makeSkewCorrectionVectors() const
 
     if (debug)
     {
-        Info<< "surfaceInterpolation::makeSkewCorrectionVectors() : "
+        Info<< "surfaceInterpolation::calcSkewCorrectionVectors() : "
             << "Finished constructing skew correction vectors"
             << endl;
     }
 }
 
 
-bool Foam::skewCorrectionVectors::skew() const
-{
-    if (skew_ == true && !skewCorrectionVectors_)
-    {
-        makeSkewCorrectionVectors();
-    }
-
-    return skew_;
-}
-
-
-const Foam::surfaceVectorField& Foam::skewCorrectionVectors::operator()() const
-{
-    if (!skew())
-    {
-        FatalErrorIn("skewCorrectionVectors::operator()()")
-            << "Cannot return skewCorrectionVectors; mesh is not skewed"
-            << abort(FatalError);
-    }
-
-    return *skewCorrectionVectors_;
-}
-
-
-// Do what is neccessary if the mesh has moved
 bool Foam::skewCorrectionVectors::movePoints()
 {
-    skew_ = true;
+    calcSkewCorrectionVectors();
-    deleteDemandDrivenData(skewCorrectionVectors_);
-
     return true;
 }
 

src/finiteVolume/interpolation/surfaceInterpolation/schemes/skewCorrected/skewCorrectionVectors.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,7 +37,7 @@ SourceFiles
 
 #include "MeshObject.H"
 #include "fvMesh.H"
-#include "surfaceFieldsFwd.H"
+#include "surfaceFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -52,18 +52,18 @@ class fvMesh;
 
 class skewCorrectionVectors
 :
-    public MeshObject<fvMesh, skewCorrectionVectors>
+    public MeshObject<fvMesh, MoveableMeshObject, skewCorrectionVectors>
 {
     // Private data
 
         //- Is mesh skew
-        mutable bool skew_;
+        bool skew_;
 
         //- Skew correction vectors
-        mutable surfaceVectorField* skewCorrectionVectors_;
+        surfaceVectorField skewCorrectionVectors_;
 
-        //- Construct skewness correction vectors
+        //- Calculate skewness correction vectors
-        void makeSkewCorrectionVectors() const;
+        void calcSkewCorrectionVectors();
 
 
 public:
@@ -83,12 +83,18 @@ public:
     // Member functions
 
         //- Return whether mesh is skew or not
-        bool skew() const;
+        bool skew() const
+        {
+            return skew_;
+        }
 
         //- Return reference to skew vectors array
-        const surfaceVectorField& operator()() const;
+        const surfaceVectorField& operator()() const
+        {
+            return skewCorrectionVectors_;
+        }
 
-        //- Delete the correction vectors when the mesh moves
+        //- Update the correction vectors when the mesh moves
         virtual bool movePoints();
 };
 

src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolation.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -612,9 +612,10 @@ void volPointInterpolation::makePatchPatchAddressing()
 
 volPointInterpolation::volPointInterpolation(const fvMesh& vm)
 :
-    MeshObject<fvMesh, volPointInterpolation>(vm)
+    MeshObject<fvMesh, Foam::UpdateableMeshObject, volPointInterpolation>(vm)
 {
-    updateMesh();
+    makeWeights();
+    makePatchPatchAddressing();
 }
 
 
@@ -626,7 +627,7 @@ volPointInterpolation::~volPointInterpolation()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void volPointInterpolation::updateMesh()
+void volPointInterpolation::updateMesh(const mapPolyMesh&)
 {
     makeWeights();
     makePatchPatchAddressing();

src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolation.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,7 +56,7 @@ class pointMesh;
 
 class volPointInterpolation
 :
-    public MeshObject<fvMesh, volPointInterpolation>
+    public MeshObject<fvMesh, UpdateableMeshObject, volPointInterpolation>
 {
     // Private data
 
@@ -174,7 +174,7 @@ public:
         // Edit
 
             //- Update mesh topology using the morph engine
-            void updateMesh();
+            void updateMesh(const mapPolyMesh&);
 
             //- Correct weighting factors for moving mesh.
             bool movePoints();

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C

@@ -120,20 +120,15 @@ void meshRefinement::testSyncBoundaryFaceList
 template<class GeoField>
 void meshRefinement::addPatchFields(fvMesh& mesh, const word& patchFieldType)
 {
-    HashTable<const GeoField*> flds
+    HashTable<GeoField*> flds
     (
         mesh.objectRegistry::lookupClass<GeoField>()
     );
 
-    forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
+    forAllIter(typename HashTable<GeoField*>, flds, iter)
     {
-        const GeoField& fld = *iter();
+        GeoField& fld = *iter();
-
+        typename GeoField::GeometricBoundaryField& bfld = fld.boundaryField();
-        typename GeoField::GeometricBoundaryField& bfld =
-            const_cast<typename GeoField::GeometricBoundaryField&>
-            (
-                fld.boundaryField()
-            );
 
         label sz = bfld.size();
         bfld.setSize(sz+1);
@@ -155,28 +150,21 @@ void meshRefinement::addPatchFields(fvMesh& mesh, const word& patchFieldType)
 template<class GeoField>
 void meshRefinement::reorderPatchFields(fvMesh& mesh, const labelList& oldToNew)
 {
-    HashTable<const GeoField*> flds
+    HashTable<GeoField*> flds
     (
         mesh.objectRegistry::lookupClass<GeoField>()
     );
 
-    forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
+    forAllIter(typename HashTable<GeoField*>, flds, iter)
     {
-        const GeoField& fld = *iter();
+        GeoField& fld = *iter();
-
+        typename GeoField::GeometricBoundaryField& bfld = fld.boundaryField();
-        typename GeoField::GeometricBoundaryField& bfld =
-            const_cast<typename GeoField::GeometricBoundaryField&>
-            (
-                fld.boundaryField()
-            );
 
         bfld.reorder(oldToNew);
     }
 }
 
 
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/meshTools/meshSearch/meshSearchFACECENTRETETSMeshObject.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -40,7 +40,12 @@ Foam::meshSearchFACECENTRETETSMeshObject::meshSearchFACECENTRETETSMeshObject
     const polyMesh& mesh
 )
 :
-    MeshObject<polyMesh, meshSearchFACECENTRETETSMeshObject>(mesh),
+    MeshObject
+    <
+        polyMesh,
+        Foam::GeometricMeshObject,
+        meshSearchFACECENTRETETSMeshObject
+    >(mesh),
     meshSearch(mesh, polyMesh::FACECENTRETETS)
 {}
 

src/meshTools/meshSearch/meshSearchFACECENTRETETSMeshObject.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -48,7 +48,12 @@ namespace Foam
 
 class meshSearchFACECENTRETETSMeshObject
 :
-    public MeshObject<polyMesh, meshSearchFACECENTRETETSMeshObject>,
+    public MeshObject
+    <
+        polyMesh,
+        GeometricMeshObject,
+        meshSearchFACECENTRETETSMeshObject
+    >,
     public meshSearch
 {
 
@@ -66,18 +71,6 @@ public:
     //- Destructor
     virtual ~meshSearchFACECENTRETETSMeshObject()
     {}
-
-//
-//    // Member functions
-//
-//        // Edit
-//
-//            //- Update mesh topology using the morph engine
-//            void updateMesh();
-//
-//            //- Correct weighting factors for moving mesh.
-//            bool movePoints();
-//
 };
 
 

src/meshTools/meshSearch/meshSearchMeshObject.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,7 +37,7 @@ namespace Foam
 
 Foam::meshSearchMeshObject::meshSearchMeshObject(const polyMesh& mesh)
 :
-    MeshObject<polyMesh, meshSearchMeshObject>(mesh),
+    MeshObject<polyMesh, Foam::GeometricMeshObject, meshSearchMeshObject>(mesh),
     meshSearch(mesh)
 {}
 

src/meshTools/meshSearch/meshSearchMeshObject.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -48,7 +48,7 @@ namespace Foam
 
 class meshSearchMeshObject
 :
-    public MeshObject<polyMesh, meshSearchMeshObject>,
+    public MeshObject<polyMesh, GeometricMeshObject, meshSearchMeshObject>,
     public meshSearch
 {
 
@@ -66,18 +66,6 @@ public:
     //- Destructor
     virtual ~meshSearchMeshObject()
     {}
-
-//
-//    // Member functions
-//
-//        // Edit
-//
-//            //- Update mesh topology using the morph engine
-//            void updateMesh();
-//
-//            //- Correct weighting factors for moving mesh.
-//            bool movePoints();
-//
 };
 
 

src/meshTools/regionSplit/regionSplit.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -692,7 +692,7 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
 
 Foam::regionSplit::regionSplit(const polyMesh& mesh)
 :
-    MeshObject<polyMesh, regionSplit>(mesh),
+    MeshObject<polyMesh, Foam::TopologicalMeshObject, regionSplit>(mesh),
     labelList(mesh.nCells(), -1)
 {
     globalNumberingPtr_ = calcRegionSplit
@@ -710,7 +710,7 @@ Foam::regionSplit::regionSplit
     const boolList& blockedFace
 )
 :
-    MeshObject<polyMesh, regionSplit>(mesh),
+    MeshObject<polyMesh, Foam::TopologicalMeshObject, regionSplit>(mesh),
     labelList(mesh.nCells(), -1)
 {
     globalNumberingPtr_ = calcRegionSplit
@@ -729,7 +729,7 @@ Foam::regionSplit::regionSplit
     const List<labelPair>& explicitConnections
 )
 :
-    MeshObject<polyMesh, regionSplit>(mesh),
+    MeshObject<polyMesh, Foam::TopologicalMeshObject, regionSplit>(mesh),
     labelList(mesh.nCells(), -1)
 {
     globalNumberingPtr_ = calcRegionSplit

src/meshTools/regionSplit/regionSplit.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -113,7 +113,7 @@ class polyMesh;
 
 class regionSplit
 :
-    public MeshObject<polyMesh, regionSplit>,
+    public MeshObject<polyMesh, TopologicalMeshObject, regionSplit>,
     public labelList
 {
     // Private data

src/parallel/decompose/scotchDecomp/scotchDecomp.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -485,6 +485,40 @@ Foam::label Foam::scotchDecomp::decomposeOneProc
             SCOTCH_archCmplt(&archdat, nProcessors_),
             "SCOTCH_archCmplt"
         );
+
+
+        //- Hack to test clustering. Note that finalDecomp is non-compact
+        //  numbers!
+        //
+        ////- Set up variable sizes architecture
+        //check
+        //(
+        //    SCOTCH_archVcmplt(&archdat),
+        //    "SCOTCH_archVcmplt"
+        //);
+        //
+        ////- stategy flags: go for quality or load balance (or leave default)
+        //SCOTCH_Num straval = 0;
+        ////straval |= SCOTCH_STRATQUALITY;
+        ////straval |= SCOTCH_STRATQUALITY;
+        //
+        ////- Number of cells per agglomeration
+        ////SCOTCH_Num agglomSize = SCOTCH_archSize(&archdat);
+        //SCOTCH_Num agglomSize = 3;
+        //
+        ////- Build strategy for agglomeration
+        //check
+        //(
+        //    SCOTCH_stratGraphClusterBuild
+        //    (
+        //        &stradat,   // strategy to build
+        //        straval,    // strategy flags
+        //        agglomSize, // cells per cluster
+        //        1.0,        // weight?
+        //        0.01        // max load imbalance
+        //    ),
+        //    "SCOTCH_stratGraphClusterBuild"
+        //);
     }
 
 

src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMeshTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -88,22 +88,20 @@ void Foam::distributedTriSurfaceMesh::distributeFields
 {
     typedef DimensionedField<Type, triSurfaceGeoMesh> DimensionedSurfField;
 
-    HashTable<const DimensionedSurfField*> fields
+    HashTable<DimensionedSurfField*> fields
     (
-        objectRegistry::lookupClass
+        objectRegistry::lookupClass<DimensionedSurfField>()
-        <DimensionedSurfField >()
     );
 
     for
     (
-        typename HashTable<const DimensionedSurfField*>::iterator fieldIter =
+        typename HashTable<DimensionedSurfField*>::iterator fieldIter =
             fields.begin();
         fieldIter != fields.end();
         ++fieldIter
     )
     {
-        DimensionedSurfField& field =
+        DimensionedSurfField& field = *fieldIter();
-            const_cast<DimensionedSurfField&>(*fieldIter());
 
         label oldSize = field.size();
 

src/postProcessing/functionObjects/field/fieldValues/fieldValueDelta/fieldValueDelta.C

@@ -27,6 +27,7 @@ License
 #include "ListOps.H"
 #include "Time.H"
 #include "volFields.H"
+#include "surfaceFields.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 

src/regionModels/pyrolysisModels/noPyrolysis/noPyrolysis.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -49,7 +49,7 @@ void noPyrolysis::constructThermoChemistry()
 {
     solidChemistry_.reset
     (
-        solidChemistryModel::New(regionMesh()).ptr()
+        basicSolidChemistryModel::New(regionMesh()).ptr()
     );
 
     solidThermo_.reset(&solidChemistry_->solidThermo());

src/regionModels/pyrolysisModels/noPyrolysis/noPyrolysis.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,7 +37,7 @@ SourceFiles
 
 #include "pyrolysisModel.H"
 #include "volFieldsFwd.H"
-#include "solidChemistryModel.H"
+#include "basicSolidChemistryModel.H"
 #include "radiationModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -82,7 +82,7 @@ protected:
         void constructThermoChemistry();
 
         //- Reference to the solid chemistry model
-        autoPtr<solidChemistryModel> solidChemistry_;
+        autoPtr<basicSolidChemistryModel> solidChemistry_;
 
         //- Reference to solid thermo
         autoPtr<solidReactionThermo> solidThermo_;

src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -303,7 +303,7 @@ void reactingOneDim::calculateMassTransfer()
 reactingOneDim::reactingOneDim(const word& modelType, const fvMesh& mesh)
 :
     pyrolysisModel(modelType, mesh),
-    solidChemistry_(solidChemistryModel::New(regionMesh())),
+    solidChemistry_(basicSolidChemistryModel::New(regionMesh())),
     solidThermo_(solidChemistry_->solidThermo()),
     radiation_(radiation::radiationModel::New(solidThermo_.T())),
     rho_
@@ -386,7 +386,7 @@ reactingOneDim::reactingOneDim
 )
 :
     pyrolysisModel(modelType, mesh, dict),
-    solidChemistry_(solidChemistryModel::New(regionMesh())),
+    solidChemistry_(basicSolidChemistryModel::New(regionMesh())),
     solidThermo_(solidChemistry_->solidThermo()),
     radiation_(radiation::radiationModel::New(solidThermo_.T())),
     rho_

src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,7 +36,7 @@ SourceFiles
 #define reactingOneDim_H
 
 #include "pyrolysisModel.H"
-#include "solidChemistryModel.H"
+#include "basicSolidChemistryModel.H"
 #include "radiationModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -76,7 +76,7 @@ protected:
     // Protected data
 
         //- Reference to the solid chemistry model
-        autoPtr<solidChemistryModel> solidChemistry_;
+        autoPtr<basicSolidChemistryModel> solidChemistry_;
 
         //- Reference to solid thermo
         solidReactionThermo& solidThermo_;

src/thermophysicalModels/SLGThermo/SLGThermo/SLGThermo.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,7 +37,15 @@ namespace Foam
 
 Foam::SLGThermo::SLGThermo(const fvMesh& mesh, fluidThermo& thermo)
 :
-    MeshObject<fvMesh, SLGThermo>(mesh),
+    regIOobject
+    (
+        IOobject
+        (
+            SLGThermo::typeName,
+            mesh.polyMesh::instance(),
+            mesh
+        )
+    ),
     thermo_(thermo),
     carrier_(NULL),
     liquids_(NULL),
@@ -246,4 +254,3 @@ bool Foam::SLGThermo::hasSolids() const
 
 
 // ************************************************************************* //
-

src/thermophysicalModels/SLGThermo/SLGThermo/SLGThermo.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -45,7 +45,7 @@ SourceFiles
 #ifndef SLGThermo_H
 #define SLGThermo_H
 
-#include "MeshObject.H"
+#include "regIOobject.H"
 #include "fluidThermo.H"
 #include "basicMultiComponentMixture.H"
 #include "liquidMixtureProperties.H"
@@ -62,7 +62,7 @@ namespace Foam
 
 class SLGThermo
 :
-    public MeshObject<fvMesh, SLGThermo>
+    public regIOobject
 {
     // Private data
 
@@ -145,6 +145,14 @@ public:
 
             //- Thermo database has solid components flag
             bool hasSolids() const;
+
+
+        // IO
+
+            bool writeData(Foam::Ostream&) const
+            {
+                return true;
+            }
 };
 
 

src/thermophysicalModels/solidChemistryModel/Make/files

@@ -1,7 +1,7 @@
 
-solidChemistryModel/solidChemistryModel.C
+basicSolidChemistryModel/basicSolidChemistryModel.C
-solidChemistryModel/solidChemistryModelNew.C
+basicSolidChemistryModel/basicSolidChemistryModelNew.C
-solidChemistryModel/solidChemistryModels.C
+basicSolidChemistryModel/basicSolidChemistryModels.C
 
 solidChemistrySolver/makeSolidChemistrySolvers.C
 

src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.C

@@ -0,0 +1,56 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "basicSolidChemistryModel.H"
+#include "fvMesh.H"
+#include "Time.H"
+
+/* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
+
+namespace Foam
+{
+    defineTypeNameAndDebug(basicSolidChemistryModel, 0);
+    defineRunTimeSelectionTable(basicSolidChemistryModel, fvMesh);
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::basicSolidChemistryModel::basicSolidChemistryModel
+(
+    const fvMesh& mesh
+)
+:
+    basicChemistryModel(mesh),
+    solidThermo_(solidReactionThermo::New(mesh))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::basicSolidChemistryModel::~basicSolidChemistryModel()
+{}
+
+
+// ************************************************************************* //

src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.H

@@ -0,0 +1,167 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013-2013 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::basicSolidChemistryModel
+
+Description
+    Chemistry model for solid thermodynamics
+
+SourceFiles
+    basicSolidChemistryModelI.H
+    basicSolidChemistryModel.C
+    newChemistrySolidModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef basicSolidChemistryModel_H
+#define basicSolidChemistryModel_H
+
+#include "basicChemistryModel.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+#include "solidReactionThermo.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class fvMesh;
+
+/*---------------------------------------------------------------------------*\
+                     class basicSolidChemistryModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class basicSolidChemistryModel
+:
+    public basicChemistryModel
+{
+    // Private Member Functions
+
+        //- Construct as copy (not implemented)
+        basicSolidChemistryModel(const basicSolidChemistryModel&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const basicSolidChemistryModel&);
+
+
+protected:
+
+    // Protected data
+
+        //- Solid thermo package
+        autoPtr<solidReactionThermo> solidThermo_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("basicSolidChemistryModel");
+
+
+    //- Declare run-time constructor selection tables
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        basicSolidChemistryModel,
+        fvMesh,
+        (const fvMesh& mesh),
+        (mesh)
+    );
+
+
+    // Constructors
+
+        //- Construct from mesh
+        basicSolidChemistryModel(const fvMesh& mesh);
+
+
+    //- Selector
+    static autoPtr<basicSolidChemistryModel> New(const fvMesh& mesh);
+
+
+    //- Destructor
+    virtual ~basicSolidChemistryModel();
+
+
+    // Member Functions
+
+        //- Return access to the solid thermo package
+        inline solidReactionThermo& solidThermo();
+
+        //- Return const access to the solid thermo package
+        inline const solidReactionThermo& solidThermo() const;
+
+        //- Return total gases mass source term [kg/m3/s]
+        virtual tmp<DimensionedField<scalar, volMesh> > RRg() const = 0;
+
+        //- Return total solids mass source term [kg/m3/s]
+        virtual tmp<DimensionedField<scalar, volMesh> > RRs() const = 0;
+
+        //- Return chemical source terms for solids [kg/m3/s]
+        virtual const DimensionedField<scalar, volMesh>& RRs
+        (
+            const label i
+        ) const = 0;
+
+        //- Return chemical source terms for gases [kg/m3/s]
+        virtual const DimensionedField<scalar, volMesh>& RRg
+        (
+            const label i
+        ) const = 0;
+
+        //- Return sensible enthalpy for gas i [J/Kg]
+        virtual tmp<volScalarField> gasHs
+        (
+            const volScalarField& p,
+            const volScalarField& T,
+            const label i
+        ) const = 0;
+
+        //- Return specie Table for gases
+        virtual const speciesTable& gasTable() const = 0;
+
+        //- Set reacting status of cell, cellI
+        virtual void setCellReacting(const label cellI, const bool active) = 0;
+
+        //- Calculates the reaction rates
+        virtual void calculate() = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "basicSolidChemistryModelI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModelI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,17 +23,19 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "fvMesh.H"
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-#include "extendedLeastSquaresGrad.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+inline Foam::solidReactionThermo& Foam::basicSolidChemistryModel::solidThermo()
+{
+    return solidThermo_();
+}
 
-namespace Foam
+
+inline const Foam::solidReactionThermo&
+Foam::basicSolidChemistryModel::solidThermo() const
 {
-namespace fv
+    return solidThermo_();
-{
-    makeFvGradScheme(extendedLeastSquaresGrad)
-}
 }
 
+
 // ************************************************************************* //

src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModelNew.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,11 +23,12 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "solidChemistryModel.H"
+#include "basicSolidChemistryModel.H"
 
 // * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
 
-Foam::autoPtr<Foam::solidChemistryModel> Foam::solidChemistryModel::New
+Foam::autoPtr<Foam::basicSolidChemistryModel> Foam::basicSolidChemistryModel::
+New
 (
     const fvMesh& mesh
 )
@@ -160,7 +161,7 @@ Foam::autoPtr<Foam::solidChemistryModel> Foam::solidChemistryModel::New
         FatalError<< exit(FatalError);
     }
 
-    return autoPtr<solidChemistryModel>(cstrIter()(mesh));
+    return autoPtr<basicSolidChemistryModel>(cstrIter()(mesh));
 }
 
 

src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModels.C

@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 InClass
-    Foam::solidChemistryModel
+    Foam::basicSolidChemistryModel
 
 Description
     Creates solid chemistry model instances templated on the type of
@@ -32,7 +32,8 @@ Description
 
 #include "makeSolidChemistryModel.H"
 
-#include "ODESolidChemistryModel.H"
+#include "pyrolysisChemistryModel.H"
+#include "basicSolidChemistryModel.H"
 #include "solidChemistryModel.H"
 #include "solidThermoPhysicsTypes.H"
 #include "thermoPhysicsTypes.H"
@@ -43,16 +44,18 @@ namespace Foam
 {
     makeSolidChemistryModel
     (
-        ODESolidChemistryModel,
         solidChemistryModel,
+        pyrolysisChemistryModel,
+        basicSolidChemistryModel,
         hConstSolidThermoPhysics,
         gasHThermoPhysics
     );
 
     makeSolidChemistryModel
     (
-        ODESolidChemistryModel,
         solidChemistryModel,
+        pyrolysisChemistryModel,
+        basicSolidChemistryModel,
         hExponentialSolidThermoPhysics,
         gasHThermoPhysics
     );

src/thermophysicalModels/solidChemistryModel/makeSolidChemistryModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,6 +30,7 @@ Description
 #define makeSolidChemistryModel_H
 
 #include "addToRunTimeSelectionTable.H"
+#include "solidChemistryModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -38,10 +39,19 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#define makeSolidChemistryModel(SS, Comp, SThermo, GThermo)                   \
+#define makeSolidChemistryModel(sChemistry, SS, Comp, SThermo, GThermo)       \
                                                                               \
     typedef SS<Comp, SThermo, GThermo> SS##Comp##SThermo##GThermo;            \
                                                                               \
+    typedef sChemistry<Comp, SThermo> sChemistryl##Comp##SThermo;             \
+                                                                              \
+    defineTemplateTypeNameAndDebugWithName                                    \
+    (                                                                         \
+        sChemistryl##Comp##SThermo,                                           \
+        (#sChemistry"<"#Comp"," + SThermo::typeName() + ">").c_str(),         \
+        0                                                                     \
+    );                                                                        \
+                                                                              \
     defineTemplateTypeNameAndDebugWithName                                    \
     (                                                                         \
         SS##Comp##SThermo##GThermo,                                           \

src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,73 +23,33 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "ODESolidChemistryModel.H"
+#include "pyrolysisChemistryModel.H"
-#include "reactingMixture.H"
 #include "solidReaction.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CompType, class SolidThermo, class GasThermo>
-Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::
+Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::
-ODESolidChemistryModel
+pyrolysisChemistryModel
 (
     const fvMesh& mesh
 )
 :
-    CompType(mesh),
+    solidChemistryModel<CompType, SolidThermo>(mesh),
-    ODE(),
+    pyrolisisGases_(this->reactions_[0].gasSpecies()),
-    Ys_(this->solidThermo().composition().Y()),
-    reactions_
-    (
-        dynamic_cast<const reactingMixture<SolidThermo>& >
-        (
-            this->solidThermo()
-        )
-    ),
-    pyrolisisGases_(reactions_[0].gasSpecies()),
-    solidThermo_
-    (
-        dynamic_cast<const reactingMixture<SolidThermo>& >
-        (
-            this->solidThermo()
-        ).speciesData()
-    ),
     gasThermo_(pyrolisisGases_.size()),
     nGases_(pyrolisisGases_.size()),
-    nSpecie_(Ys_.size() + nGases_),
+    nSpecie_(this->Ys_.size() + nGases_),
-    nSolids_(Ys_.size()),
-    nReaction_(reactions_.size()),
-    RRs_(nSolids_),
     RRg_(nGases_),
-    Ys0_(nSolids_),
+    Ys0_(this->nSolids_),
-    cellCounter_(0),
+    cellCounter_(0)
-    reactingCells_(mesh.nCells(), true)
 {
     // create the fields for the chemistry sources
-    forAll(RRs_, fieldI)
+    forAll(this->RRs_, fieldI)
     {
-        RRs_.set
-        (
-            fieldI,
-            new DimensionedField<scalar, volMesh>
-            (
-                IOobject
-                (
-                    "RRs." + Ys_[fieldI].name(),
-                    mesh.time().timeName(),
-                    mesh,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
-                mesh,
-                dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
-            )
-        );
-
-
         IOobject header
         (
-            Ys_[fieldI].name() + "0",
+            this->Ys_[fieldI].name() + "0",
             mesh.time().timeName(),
             mesh,
             IOobject::NO_READ
@@ -105,7 +65,7 @@ ODESolidChemistryModel
                 (
                     IOobject
                     (
-                        Ys_[fieldI].name() + "0",
+                        this->Ys_[fieldI].name() + "0",
                         mesh.time().timeName(),
                         mesh,
                         IOobject::MUST_READ,
@@ -137,7 +97,7 @@ ODESolidChemistryModel
                 (
                     IOobject
                     (
-                        Ys_[fieldI].name() + "0",
+                        this->Ys_[fieldI].name() + "0",
                         mesh.time().timeName(),
                         mesh,
                         IOobject::NO_READ,
@@ -150,7 +110,7 @@ ODESolidChemistryModel
             // Calculate inital values of Ysi0 = rho*delta*Yi
             Ys0_[fieldI].internalField() =
                 this->solidThermo().rho()
-               *max(Ys_[fieldI], scalar(0.001))*mesh.V();
+               *max(this->Ys_[fieldI], scalar(0.001))*mesh.V();
         }
     }
 
@@ -190,23 +150,24 @@ ODESolidChemistryModel
         );
     }
 
-    Info<< "solidChemistryModel: Number of solids = " << nSolids_ << endl;
+    Info<< "pyrolysisChemistryModel: " << nl;
-
+    Info<< indent << "Number of solids = " << this->nSolids_ << nl;
-    Info<< "solidChemistryModel: Number of gases = " << nGases_ << endl;
+    Info<< indent << "Number of gases = " << nGases_ << nl;
-
+    forAll(this->reactions_, i)
-    forAll(reactions_, i)
     {
-        Info<< indent << reactions_[i] << nl;
+        Info<< dynamic_cast<const solidReaction<SolidThermo>& >
+        (
+            this->reactions_[i]
+        ) << nl;
     }
-
 }
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 template<class CompType, class SolidThermo, class GasThermo>
-Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::
+Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::
-~ODESolidChemistryModel()
+~pyrolysisChemistryModel()
 {}
 
 
@@ -214,7 +175,7 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::
 
 template<class CompType, class SolidThermo, class GasThermo>
 Foam::scalarField Foam::
-ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::omega
+pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::omega
 (
     const scalarField& c,
     const scalar T,
@@ -229,9 +190,9 @@ ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::omega
 
     scalarField om(nEqns(), 0.0);
 
-    forAll(reactions_, i)
+    forAll(this->reactions_, i)
     {
-        const Reaction<SolidThermo>& R = reactions_[i];
+        const Reaction<SolidThermo>& R = this->reactions_[i];
 
         scalar omegai = omega
         (
@@ -242,13 +203,13 @@ ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::omega
         {
             label si = R.lhs()[s].index;
             om[si] -= omegai;
-            rhoL = solidThermo_[si].rho(p, T);
+            rhoL = this->solidThermo_[si].rho(p, T);
         }
         scalar sr = 0.0;
         forAll(R.rhs(), s)
         {
             label si = R.rhs()[s].index;
-            scalar rhoR = solidThermo_[si].rho(p, T);
+            scalar rhoR = this->solidThermo_[si].rho(p, T);
             sr = rhoR/rhoL;
             om[si] += sr*omegai;
 
@@ -260,7 +221,7 @@ ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::omega
         forAll(R.grhs(), g)
         {
             label gi = R.grhs()[g].index;
-            om[gi + nSolids_] += (1.0 - sr)*omegai;
+            om[gi + this->nSolids_] += (1.0 - sr)*omegai;
         }
     }
 
@@ -270,7 +231,7 @@ ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::omega
 
 template<class CompType, class SolidThermo, class GasThermo>
 Foam::scalar
-Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::omega
+Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::omega
 (
     const Reaction<SolidThermo>& R,
     const scalarField& c,
@@ -314,7 +275,7 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::omega
 
 
 template<class CompType, class SolidThermo, class GasThermo>
-void Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::
+void Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::
 derivatives
 (
     const scalar time,
@@ -331,19 +292,19 @@ derivatives
 
     //Total mass concentration
     scalar cTot = 0.0;
-    for (label i=0; i<nSolids_; i++)
+    for (label i=0; i<this->nSolids_; i++)
     {
         cTot += c[i];
     }
 
     scalar newCp = 0.0;
     scalar newhi = 0.0;
-    for (label i=0; i<nSolids_; i++)
+    for (label i=0; i<this->nSolids_; i++)
     {
         scalar dYidt = dcdt[i]/cTot;
         scalar Yi = c[i]/cTot;
-        newCp += Yi*solidThermo_[i].Cp(p, T);
+        newCp += Yi*this->solidThermo_[i].Cp(p, T);
-        newhi -= dYidt*solidThermo_[i].Hc();
+        newhi -= dYidt*this->solidThermo_[i].Hc();
     }
 
     scalar dTdt = newhi/newCp;
@@ -356,7 +317,8 @@ derivatives
 
 
 template<class CompType, class SolidThermo, class GasThermo>
-void Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::jacobian
+void Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::
+jacobian
 (
     const scalar t,
     const scalarField& c,
@@ -369,7 +331,7 @@ void Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::jacobian
 
     scalarField c2(nSpecie_, 0.0);
 
-    for (label i=0; i<nSolids_; i++)
+    for (label i=0; i<this->nSolids_; i++)
     {
         c2[i] = max(c[i], 0.0);
     }
@@ -385,9 +347,9 @@ void Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::jacobian
     // length of the first argument must be nSolids
     dcdt = omega(c2, T, p);
 
-    for (label ri=0; ri<reactions_.size(); ri++)
+    for (label ri=0; ri<this->reactions_.size(); ri++)
     {
-        const Reaction<SolidThermo>& R = reactions_[ri];
+        const Reaction<SolidThermo>& R = this->reactions_[ri];
 
         scalar kf0 = R.kf(p, T, c2);
 
@@ -451,96 +413,9 @@ void Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::jacobian
 }
 
 
-template<class CompType, class SolidThermo, class GasThermo>
-Foam::tmp<Foam::volScalarField>
-Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::tc() const
-{
-    notImplemented
-    (
-        "ODESolidChemistryModel::tc()"
-    );
-
-    return volScalarField::null();
-}
-
-
-template<class CompType, class SolidThermo, class GasThermo>
-Foam::tmp<Foam::volScalarField>
-Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::Sh() const
-{
-    tmp<volScalarField> tSh
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "Sh",
-                this->mesh_.time().timeName(),
-                this->mesh_,
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE,
-                false
-            ),
-            this->mesh_,
-            dimensionedScalar("zero", dimEnergy/dimTime/dimVolume, 0.0),
-            zeroGradientFvPatchScalarField::typeName
-        )
-    );
-
-    if (this->chemistry_)
-    {
-        scalarField& Sh = tSh();
-
-        forAll(Ys_, i)
-        {
-            forAll(Sh, cellI)
-            {
-                scalar hf = solidThermo_[i].Hc();
-                Sh[cellI] -= hf*RRs_[i][cellI];
-            }
-        }
-    }
-
-    return tSh;
-}
-
-
-template<class CompType, class SolidThermo, class GasThermo>
-Foam::tmp<Foam::volScalarField>
-Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::dQ() const
-{
-    tmp<volScalarField> tdQ
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "dQ",
-                this->mesh_.time().timeName(),
-                this->mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            this->mesh_,
-            dimensionedScalar("dQ", dimEnergy/dimTime, 0.0),
-            zeroGradientFvPatchScalarField::typeName
-        )
-    );
-
-    if (this->chemistry_)
-    {
-        volScalarField& dQ = tdQ();
-        dQ.dimensionedInternalField() = this->mesh_.V()*Sh()();
-    }
-
-    return tdQ;
-}
-
-
 template<class CompType, class SolidThermo, class GasThermo>
 Foam::label Foam::
-ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::nEqns() const
+pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::nEqns() const
 {
     // nEqns = number of solids + gases + temperature + pressure
     return (nSpecie_ + 2);
@@ -548,7 +423,7 @@ ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::nEqns() const
 
 
 template<class CompType, class SolidThermo, class GasThermo>
-void Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::
+void Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::
 calculate()
 {
     if (!this->chemistry_)
@@ -570,10 +445,11 @@ calculate()
         this->solidThermo().rho()
     );
 
-    forAll(RRs_, i)
+    forAll(this->RRs_, i)
     {
-        RRs_[i].field() = 0.0;
+        this->RRs_[i].field() = 0.0;
     }
+
     forAll(RRg_, i)
     {
         RRg_[i].field() = 0.0;
@@ -585,28 +461,28 @@ calculate()
 
         const scalar delta = this->mesh().V()[celli];
 
-        if (reactingCells_[celli])
+        if (this->reactingCells_[celli])
         {
             scalar rhoi = rho[celli];
             scalar Ti = this->solidThermo().T()[celli];
             scalar pi = this->solidThermo().p()[celli];
 
             scalarField c(nSpecie_, 0.0);
-            for (label i=0; i<nSolids_; i++)
+            for (label i=0; i<this->nSolids_; i++)
             {
-                c[i] = rhoi*Ys_[i][celli]*delta;
+                c[i] = rhoi*this->Ys_[i][celli]*delta;
             }
 
             const scalarField dcdt = omega(c, Ti, pi, true);
 
-            forAll(RRs_, i)
+            forAll(this->RRs_, i)
             {
-                RRs_[i][celli] = dcdt[i]/delta;
+                this->RRs_[i][celli] = dcdt[i]/delta;
             }
 
             forAll(RRg_, i)
             {
-                RRg_[i][celli] = dcdt[nSolids_ + i]/delta;
+                RRg_[i][celli] = dcdt[this->nSolids_ + i]/delta;
             }
         }
     }
@@ -615,7 +491,7 @@ calculate()
 
 template<class CompType, class SolidThermo, class GasThermo>
 Foam::scalar
-Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::solve
+Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::solve
 (
     const scalar t0,
     const scalar deltaT
@@ -642,9 +518,9 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::solve
         this->solidThermo().rho()
     );
 
-    forAll(RRs_, i)
+    forAll(this->RRs_, i)
     {
-        RRs_[i].field() = 0.0;
+        this->RRs_[i].field() = 0.0;
     }
     forAll(RRg_, i)
     {
@@ -654,7 +530,7 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::solve
 
     forAll(rho, celli)
     {
-        if (reactingCells_[celli])
+        if (this->reactingCells_[celli])
         {
             cellCounter_ = celli;
 
@@ -668,9 +544,9 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::solve
 
             scalar delta = this->mesh().V()[celli];
 
-            for (label i=0; i<nSolids_; i++)
+            for (label i=0; i<this->nSolids_; i++)
             {
-                c[i] = rhoi*Ys_[i][celli]*delta;
+                c[i] = rhoi*this->Ys_[i][celli]*delta;
             }
 
             c0 = c;
@@ -690,7 +566,7 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::solve
                 scalar cTot = 0.0;
 
                 //Total mass concentration
-                for (label i=0; i<nSolids_; i++)
+                for (label i=0; i<this->nSolids_; i++)
                 {
                     cTot += c[i];
                 }
@@ -700,13 +576,13 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::solve
                 scalar invRho = 0.0;
                 scalarList dcdt = (c - c0)/dt;
 
-                for (label i=0; i<nSolids_; i++)
+                for (label i=0; i<this->nSolids_; i++)
                 {
                     scalar dYi = dcdt[i]/cTot;
                     scalar Yi = c[i]/cTot;
-                    newCp += Yi*solidThermo_[i].Cp(pi, Ti);
+                    newCp += Yi*this->solidThermo_[i].Cp(pi, Ti);
-                    newhi -= dYi*solidThermo_[i].Hc();
+                    newhi -= dYi*this->solidThermo_[i].Hc();
-                    invRho += Yi/solidThermo_[i].rho(pi, Ti);
+                    invRho += Yi/this->solidThermo_[i].rho(pi, Ti);
                 }
 
                 scalar dTi = (newhi/newCp)*dt;
@@ -722,14 +598,14 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::solve
             deltaTMin = min(tauC, deltaTMin);
             dc = c - c0;
 
-            forAll(RRs_, i)
+            forAll(this->RRs_, i)
             {
-                RRs_[i][celli] = dc[i]/(deltaT*delta);
+                this->RRs_[i][celli] = dc[i]/(deltaT*delta);
             }
 
             forAll(RRg_, i)
             {
-                RRg_[i][celli] = dc[nSolids_ + i]/(deltaT*delta);
+                RRg_[i][celli] = dc[this->nSolids_ + i]/(deltaT*delta);
             }
 
             // Update Ys0_
@@ -746,7 +622,7 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::solve
 
 template<class CompType, class SolidThermo,class GasThermo>
 Foam::tmp<Foam::volScalarField>
-Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::gasHs
+Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::gasHs
 (
     const volScalarField& p,
     const volScalarField& T,
@@ -785,9 +661,9 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::gasHs
 }
 
 
-template<class CompType, class SolidThermo,class GasThermo>
+template<class CompType, class SolidThermo, class GasThermo>
 Foam::scalar
-Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::solve
+Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::solve
 (
     scalarField &c,
     const scalar T,
@@ -798,7 +674,7 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::solve
 {
     notImplemented
     (
-        "ODESolidChemistryModel::solve"
+        "pyrolysisChemistryModel::solve"
         "("
             "scalarField&, "
             "const scalar, "
@@ -809,14 +685,4 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::solve
     );
     return (0);
 }
-
-
-template<class CompType, class SolidThermo,class GasThermo>
-void Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::
-setCellReacting(const label cellI, const bool active)
-{
-    reactingCells_[cellI] = active;
-}
-
-
 // ************************************************************************* //

src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,26 +22,24 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::ODESolidChemistryModel
+    Foam::pyrolysisChemistryModel
 
 Description
-    Extends base chemistry model by adding a thermo package, and ODE functions.
+    Pyrolysis chemistry model. It includes gas phase in the solid
-    Introduces chemistry equation system and evaluation of chemical source
+    reaction.
-    terms.
 
 SourceFiles
-    ODESolidChemistryModelI.H
+    pyrolysisChemistryModelI.H
-    ODESolidChemistryModel.C
+    pyrolysisChemistryModel.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef ODESolidChemistryModel_H
+#ifndef pyrolysisChemistryModel_H
-#define ODESolidChemistryModel_H
+#define pyrolysisChemistryModel_H
 
-#include "Reaction.H"
-#include "ODE.H"
 #include "volFieldsFwd.H"
 #include "DimensionedField.H"
+#include "solidChemistryModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -52,35 +50,25 @@ namespace Foam
 class fvMesh;
 
 /*---------------------------------------------------------------------------*\
-                   Class ODESolidChemistryModel Declaration
+                   Class pyrolysisChemistryModel Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class CompType, class SolidThermo, class GasThermo>
-class ODESolidChemistryModel
+class pyrolysisChemistryModel
 :
-    public CompType,
+    public solidChemistryModel<CompType, SolidThermo>
-    public ODE
 {
     // Private Member Functions
 
         //- Disallow default bitwise assignment
-        void operator=(const ODESolidChemistryModel&);
+        void operator=(const pyrolysisChemistryModel&);
 
 
 protected:
 
-        //- Reference to solid mass fractions
-        PtrList<volScalarField>& Ys_;
-
-        //- Reactions
-        const PtrList<Reaction<SolidThermo> >& reactions_;
-
         //- List of gas species present in reaction system
         speciesTable pyrolisisGases_;
 
-        //- Thermodynamic data of solids
-        const PtrList<SolidThermo>& solidThermo_;
-
         //- Thermodynamic data of gases
         PtrList<GasThermo> gasThermo_;
 
@@ -90,24 +78,12 @@ protected:
         //- Number of components being solved by ODE
         label nSpecie_;
 
-        //- Number of solid components
-        label nSolids_;
-
-        //- Number of solid reactions
-        label nReaction_;
-
-        //- List of reaction rate per solid [kg/m3/s]
-        PtrList<DimensionedField<scalar, volMesh> > RRs_;
-
         //- List of reaction rate per gas [kg/m3/s]
         PtrList<DimensionedField<scalar, volMesh> > RRg_;
 
 
     // Protected Member Functions
 
-        //- Write access to source terms for solids
-        inline PtrList<DimensionedField<scalar, volMesh> >& RRs();
-
         //- Write access to source terms for gases
         inline PtrList<DimensionedField<scalar, volMesh> >& RRg();
 
@@ -120,37 +96,25 @@ private:
         //- Cell counter
         label cellCounter_;
 
-        //- List of active reacting cells
-        List<bool> reactingCells_;
-
-
-        // Private members
-
-            //- Set reacting status of cell, cellI
-            void setCellReacting(const label cellI, const bool active);
-
 
 public:
 
     //- Runtime type information
-    TypeName("ODESolidChemistryModel");
+    TypeName("pyrolysis");
 
 
     // Constructors
 
         //- Construct from mesh
-        ODESolidChemistryModel(const fvMesh& mesh);
+        pyrolysisChemistryModel(const fvMesh& mesh);
 
 
     //- Destructor
-    virtual ~ODESolidChemistryModel();
+    virtual ~pyrolysisChemistryModel();
 
 
     // Member Functions
 
-        //- The reactions
-        inline const PtrList<Reaction<SolidThermo> >& reactions() const;
-
         //- Thermodynamic data of gases
         inline const PtrList<GasThermo>& gasThermo() const;
 
@@ -163,9 +127,6 @@ public:
         //- The number of solids
         inline label nGases() const;
 
-        //- The number of reactions
-        inline label nReaction() const;
-
 
         //- dc/dt = omega, rate of change in concentration, for each species
         virtual scalarField omega
@@ -198,12 +159,6 @@ public:
 
         // Chemistry model functions
 
-            //- Return const access to the chemical source terms for solids
-            inline const DimensionedField<scalar, volMesh>& RRs
-            (
-                const label i
-            ) const;
-
             //- Return const access to the chemical source terms for gases
             inline const DimensionedField<scalar, volMesh>& RRg
             (
@@ -213,15 +168,6 @@ public:
             //- Return total gas source term
             inline tmp<DimensionedField<scalar, volMesh> > RRg() const;
 
-            //- Return total solid source term
-            inline tmp<DimensionedField<scalar, volMesh> > RRs() const;
-
-            //- Return const access to the total source terms
-            inline const DimensionedField<scalar, volMesh>& RR
-            (
-                const label i
-            ) const;
-
             //- Return sensible enthalpy for gas i [J/Kg]
             virtual tmp<volScalarField> gasHs
             (
@@ -232,16 +178,7 @@ public:
 
             //- Solve the reaction system for the given start time and time
             //  step and return the characteristic time
-            virtual scalar solve(const scalar t0, const scalar deltaT);
+            virtual scalar solve(const scalar t0, const scalar deltaT) ;
-
-            //- Return the chemical time scale
-            virtual tmp<volScalarField> tc() const;
-
-            //- Return source for enthalpy equation [kg/m/s3]
-            virtual tmp<volScalarField> Sh() const;
-
-            //- Return the heat release, i.e. enthalpy/sec [m2/s3]
-            virtual tmp<volScalarField> dQ() const;
 
 
         // ODE functions (overriding abstract functions in ODE.H)
@@ -281,12 +218,12 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#   include "ODESolidChemistryModelI.H"
+#   include "pyrolysisChemistryModelI.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-#   include "ODESolidChemistryModel.C"
+#   include "pyrolysisChemistryModel.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModelI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,37 +24,21 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "volFields.H"
-#include "zeroGradientFvPatchFields.H"
+
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CompType, class SolidThermo, class GasThermo>
 inline Foam::PtrList<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >&
-Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RRs()
+Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::RRg()
-{
-    return RRs_;
-}
-
-
-template<class CompType, class SolidThermo, class GasThermo>
-inline Foam::PtrList<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >&
-Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RRg()
 {
     return RRg_;
 }
 
 
-template<class CompType, class SolidThermo, class GasThermo>
-inline const Foam::PtrList<Foam::Reaction<SolidThermo> >&
-Foam::ODESolidChemistryModel<CompType, SolidThermo,GasThermo>::reactions() const
-{
-    return reactions_;
-}
-
-
 template<class CompType, class SolidThermo, class GasThermo>
 inline const Foam::PtrList<GasThermo>&
-Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::
+Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::
 gasThermo() const
 {
     return gasThermo_;
@@ -63,7 +47,8 @@ gasThermo() const
 
 template<class CompType, class SolidThermo, class GasThermo>
 inline const Foam::speciesTable&
-Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::gasTable() const
+Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::
+gasTable() const
 {
     return pyrolisisGases_;
 }
@@ -71,35 +56,16 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::gasTable() const
 
 template<class CompType, class SolidThermo, class GasThermo>
 inline Foam::label
-Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::nSpecie() const
+Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::
+nSpecie() const
 {
     return nSpecie_;
 }
 
 
-template<class CompType, class SolidThermo, class GasThermo>
-inline Foam::label
-Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::
-nReaction() const
-{
-    return nReaction_;
-}
-
-
 template<class CompType, class SolidThermo, class GasThermo>
 inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
-Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RRs
+Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::RRg
-(
-    const label i
-) const
-{
-    return RRs_[i];
-}
-
-
-template<class CompType, class SolidThermo, class GasThermo>
-inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
-Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RRg
 (
     const label i
 ) const
@@ -110,7 +76,8 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RRg
 
 template<class CompType, class SolidThermo, class GasThermo>
 inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
-Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RRg() const
+Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::
+RRg() const
 {
     tmp<DimensionedField<scalar, volMesh> > tRRg
     (
@@ -139,51 +106,4 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RRg() const
     }
     return tRRg;
 }
-
-
-template<class CompType, class SolidThermo, class GasThermo>
-inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
-Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RRs() const
-{
-    tmp<DimensionedField<scalar, volMesh> > tRRs
-    (
-        new DimensionedField<scalar, volMesh>
-        (
-            IOobject
-            (
-                "RRs",
-                this->time().timeName(),
-                this->mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            this->mesh(),
-            dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
-        )
-    );
-
-    if (this->chemistry_)
-    {
-        DimensionedField<scalar, volMesh>& RRs = tRRs();
-        for (label i=0; i < nSolids_; i++)
-        {
-            RRs += RRs_[i];
-        }
-    }
-    return tRRs;
-}
-
-
-template<class CompType, class SolidThermo, class GasThermo>
-inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
-Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RR
-(
-    const label i
-) const
-{
-    notImplemented("ODESolidChemistryModel::RR(const label)");
-    return (DimensionedField<scalar, volMesh>::null());
-}
-
-
 // ************************************************************************* //

src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,33 +24,193 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "solidChemistryModel.H"
-#include "fvMesh.H"
+#include "reactingMixture.H"
-#include "Time.H"
+#include "zeroGradientFvPatchFields.H"
-
-/* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
-
-namespace Foam
-{
-    defineTypeNameAndDebug(solidChemistryModel, 0);
-    defineRunTimeSelectionTable(solidChemistryModel, fvMesh);
-}
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::solidChemistryModel::solidChemistryModel
+template<class CompType, class SolidThermo>
+Foam::solidChemistryModel<CompType, SolidThermo>::
+solidChemistryModel
 (
     const fvMesh& mesh
 )
 :
-    basicChemistryModel(mesh),
+    CompType(mesh),
-    solidThermo_(solidReactionThermo::New(mesh))
+    ODE(),
-{}
+    Ys_(this->solidThermo().composition().Y()),
+    reactions_
+    (
+        dynamic_cast<const reactingMixture<SolidThermo>& >
+        (
+            this->solidThermo()
+        )
+    ),
+    solidThermo_
+    (
+        dynamic_cast<const reactingMixture<SolidThermo>& >
+        (
+            this->solidThermo()
+        ).speciesData()
+    ),
+    nSolids_(Ys_.size()),
+    nReaction_(reactions_.size()),
+    RRs_(nSolids_),
+    reactingCells_(mesh.nCells(), true)
+{
+    // create the fields for the chemistry sources
+    forAll(RRs_, fieldI)
+    {
+        RRs_.set
+        (
+            fieldI,
+            new DimensionedField<scalar, volMesh>
+            (
+                IOobject
+                (
+                    "RRs." + Ys_[fieldI].name(),
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh,
+                dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
+            )
+        );
+   }
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-Foam::solidChemistryModel::~solidChemistryModel()
+template<class CompType, class SolidThermo>
+Foam::solidChemistryModel<CompType, SolidThermo>::
+~solidChemistryModel()
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CompType, class SolidThermo>
+Foam::tmp<Foam::volScalarField>
+Foam::solidChemistryModel<CompType, SolidThermo>::tc() const
+{
+    notImplemented
+    (
+        "solidChemistryModel::tc()"
+    );
+    return volScalarField::null();
+}
+
+
+template<class CompType, class SolidThermo>
+Foam::tmp<Foam::volScalarField>
+Foam::solidChemistryModel<CompType, SolidThermo>::Sh() const
+{
+    tmp<volScalarField> tSh
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Sh",
+                this->mesh_.time().timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            this->mesh_,
+            dimensionedScalar("zero", dimEnergy/dimTime/dimVolume, 0.0),
+            zeroGradientFvPatchScalarField::typeName
+        )
+    );
+
+    if (this->chemistry_)
+    {
+        scalarField& Sh = tSh();
+
+        forAll(Ys_, i)
+        {
+            forAll(Sh, cellI)
+            {
+                scalar hf = solidThermo_[i].Hc();
+                Sh[cellI] -= hf*RRs_[i][cellI];
+            }
+        }
+    }
+
+    return tSh;
+}
+
+
+template<class CompType, class SolidThermo>
+Foam::tmp<Foam::volScalarField>
+Foam::solidChemistryModel<CompType, SolidThermo>::dQ() const
+{
+    tmp<volScalarField> tdQ
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "dQ",
+                this->mesh_.time().timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            this->mesh_,
+            dimensionedScalar("dQ", dimEnergy/dimTime, 0.0),
+            zeroGradientFvPatchScalarField::typeName
+        )
+    );
+
+    if (this->chemistry_)
+    {
+        volScalarField& dQ = tdQ();
+        dQ.dimensionedInternalField() = this->mesh_.V()*Sh()();
+    }
+
+    return tdQ;
+}
+
+
+template<class CompType, class SolidThermo>
+Foam::scalar Foam::solidChemistryModel<CompType, SolidThermo>::solve
+(
+    scalarField &c,
+    const scalar T,
+    const scalar p,
+    const scalar t0,
+    const scalar dt
+) const
+{
+    notImplemented
+    (
+        "solidChemistryModel::solve"
+        "("
+            "scalarField&, "
+            "const scalar, "
+            "const scalar, "
+            "const scalar, "
+            "const scalar"
+        ") const"
+    );
+    return (0);
+}
+
+
+template<class CompType, class SolidThermo>
+void Foam::solidChemistryModel<CompType, SolidThermo>::setCellReacting
+(
+    const label cellI,
+    const bool active
+)
+{
+    reactingCells_[cellI] = active;
+}
+
 // ************************************************************************* //

src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,23 +25,24 @@ Class
     Foam::solidChemistryModel
 
 Description
-    Chemistry model for solid thermodynamics
+    Extends base solid chemistry model by adding a thermo package, and ODE
+    functions.
+    Introduces chemistry equation system and evaluation of chemical source
+    terms.
 
 SourceFiles
     solidChemistryModelI.H
     solidChemistryModel.C
-    newChemistrySolidModel.C
 
 \*---------------------------------------------------------------------------*/
 
 #ifndef solidChemistryModel_H
 #define solidChemistryModel_H
 
-#include "basicChemistryModel.H"
+#include "Reaction.H"
-#include "autoPtr.H"
+#include "ODE.H"
-#include "runTimeSelectionTables.H"
+#include "volFieldsFwd.H"
-#include "solidReactionThermo.H"
+#include "DimensionedField.H"
-
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -52,45 +53,57 @@ namespace Foam
 class fvMesh;
 
 /*---------------------------------------------------------------------------*\
-                     class solidChemistryModel Declaration
+                   Class solidChemistryModel Declaration
 \*---------------------------------------------------------------------------*/
 
+template<class CompType, class SolidThermo>
 class solidChemistryModel
 :
-    public basicChemistryModel
+    public CompType,
+    public ODE
 {
     // Private Member Functions
 
-        //- Construct as copy (not implemented)
-        solidChemistryModel(const solidChemistryModel&);
-
         //- Disallow default bitwise assignment
         void operator=(const solidChemistryModel&);
 
 
 protected:
 
-    // Protected data
+        //- Reference to solid mass fractions
+        PtrList<volScalarField>& Ys_;
 
-        //- Solid thermo package
+        //- Reactions
-        autoPtr<solidReactionThermo> solidThermo_;
+        const PtrList<Reaction<SolidThermo> >& reactions_;
+
+        //- Thermodynamic data of solids
+        const PtrList<SolidThermo>& solidThermo_;
+
+        //- Number of solid components
+        label nSolids_;
+
+        //- Number of solid reactions
+        label nReaction_;
+
+        //- List of reaction rate per solid [kg/m3/s]
+        PtrList<DimensionedField<scalar, volMesh> > RRs_;
+
+        //- List of active reacting cells
+        List<bool> reactingCells_;
+
+    // Protected Member Functions
+
+        //- Write access to source terms for solids
+        inline PtrList<DimensionedField<scalar, volMesh> >& RRs();
+
+        //- Set reacting status of cell, cellI
+        void setCellReacting(const label cellI, const bool active);
 
 
 public:
 
     //- Runtime type information
-    TypeName("solid");
+    TypeName("solidChemistryModel");
-
-
-    //- Declare run-time constructor selection tables
-    declareRunTimeSelectionTable
-    (
-        autoPtr,
-        solidChemistryModel,
-        fvMesh,
-        (const fvMesh& mesh),
-        (mesh)
-    );
 
 
     // Constructors
@@ -99,56 +112,108 @@ public:
         solidChemistryModel(const fvMesh& mesh);
 
 
-    //- Selector
-    static autoPtr<solidChemistryModel> New(const fvMesh& mesh);
-
-
     //- Destructor
     virtual ~solidChemistryModel();
 
 
     // Member Functions
 
-        //- Return access to the solid thermo package
+        //- The reactions
-        inline solidReactionThermo& solidThermo();
+        inline const PtrList<Reaction<SolidThermo> >& reactions() const;
 
-        //- Return const access to the solid thermo package
+        //- The number of reactions
-        inline const solidReactionThermo& solidThermo() const;
+        inline label nReaction() const;
 
-        //- Return total gases mass source term [kg/m3/s]
-        virtual tmp<DimensionedField<scalar, volMesh> > RRg() const = 0;
 
-        //- Return total solids mass source term [kg/m3/s]
+        //- dc/dt = omega, rate of change in concentration, for each species
-        virtual tmp<DimensionedField<scalar, volMesh> > RRs() const = 0;
+        virtual scalarField omega
-
-        //- Return chemical source terms for solids [kg/m3/s]
-        virtual const DimensionedField<scalar, volMesh>& RRs
         (
-            const label i
+            const scalarField& c,
+            const scalar T,
+            const scalar p,
+            const bool updateC0 = false
         ) const = 0;
 
-        //- Return chemical source terms for gases [kg/m3/s]
+        //- Return the reaction rate for reaction r and the reference
-        virtual const DimensionedField<scalar, volMesh>& RRg
+        //  species and charateristic times
+        virtual scalar omega
         (
-            const label i
+            const Reaction<SolidThermo>& r,
+            const scalarField& c,
+            const scalar T,
+            const scalar p,
+            scalar& pf,
+            scalar& cf,
+            label& lRef,
+            scalar& pr,
+            scalar& cr,
+            label& rRef
         ) const = 0;
 
-        //- Return sensible enthalpy for gas i [J/Kg]
-        virtual tmp<volScalarField> gasHs
-        (
-            const volScalarField& p,
-            const volScalarField& T,
-            const label i
-        ) const = 0;
-
-        //- Return specie Table for gases
-        virtual const speciesTable& gasTable() const = 0;
-
-        //- Set reacting status of cell, cellI
-        virtual void setCellReacting(const label cellI, const bool active) = 0;
-
         //- Calculates the reaction rates
         virtual void calculate() = 0;
+
+
+        // Chemistry model functions
+
+            //- Return const access to the chemical source terms for solids
+            inline const DimensionedField<scalar, volMesh>& RRs
+            (
+                const label i
+            ) const;
+
+            //- Return total solid source term
+            inline tmp<DimensionedField<scalar, volMesh> > RRs() const;
+
+            //- Return const access to the total source terms
+            inline const DimensionedField<scalar, volMesh>& RR
+            (
+                const label i
+            ) const;
+
+
+            //- Solve the reaction system for the given start time and time
+            //  step and return the characteristic time
+            virtual scalar solve(const scalar t0, const scalar deltaT) = 0;
+
+            //- Return the chemical time scale
+            virtual tmp<volScalarField> tc() const;
+
+            //- Return source for enthalpy equation [kg/m/s3]
+            virtual tmp<volScalarField> Sh() const;
+
+            //- Return the heat release, i.e. enthalpy/sec [m2/s3]
+            virtual tmp<volScalarField> dQ() const;
+
+
+        // ODE functions (overriding abstract functions in ODE.H)
+
+            //- Number of ODE's to solve
+            virtual label nEqns() const = 0;
+
+            virtual void derivatives
+            (
+                const scalar t,
+                const scalarField& c,
+                scalarField& dcdt
+            ) const = 0;
+
+            virtual void jacobian
+            (
+                const scalar t,
+                const scalarField& c,
+                scalarField& dcdt,
+                scalarSquareMatrix& dfdc
+            ) const = 0;
+
+            virtual scalar solve
+            (
+                scalarField &c,
+                const scalar T,
+                const scalar p,
+                const scalar t0,
+                const scalar dt
+            ) const = 0;
 };
 
 
@@ -158,7 +223,13 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#include "solidChemistryModelI.H"
+#   include "solidChemistryModelI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "solidChemistryModel.C"
+#endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModelI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,18 +23,87 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+#include "volFields.H"
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline Foam::solidReactionThermo& Foam::solidChemistryModel::solidThermo()
+template<class CompType, class SolidThermo>
+inline Foam::PtrList<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >&
+Foam::solidChemistryModel<CompType, SolidThermo>::RRs()
 {
-    return solidThermo_();
+    return RRs_;
+}
+
+template<class CompType, class SolidThermo>
+inline const Foam::PtrList<Foam::Reaction<SolidThermo> >&
+Foam::solidChemistryModel<CompType, SolidThermo>::reactions() const
+{
+    return reactions_;
 }
 
 
-inline const Foam::solidReactionThermo&
+template<class CompType, class SolidThermo>
-Foam::solidChemistryModel::solidThermo() const
+inline Foam::label
+Foam::solidChemistryModel<CompType, SolidThermo>::
+nReaction() const
 {
-    return solidThermo_();
+    return nReaction_;
+}
+
+
+template<class CompType, class SolidThermo>
+inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::solidChemistryModel<CompType, SolidThermo>::RRs
+(
+    const label i
+) const
+{
+    return RRs_[i];
+}
+
+
+template<class CompType, class SolidThermo>
+inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+Foam::solidChemistryModel<CompType, SolidThermo>::RRs() const
+{
+    tmp<DimensionedField<scalar, volMesh> > tRRs
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                "RRs",
+                this->time().timeName(),
+                this->mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            this->mesh(),
+            dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
+        )
+    );
+
+    if (this->chemistry_)
+    {
+        DimensionedField<scalar, volMesh>& RRs = tRRs();
+        for (label i=0; i < nSolids_; i++)
+        {
+            RRs += RRs_[i];
+        }
+    }
+    return tRRs;
+}
+
+
+template<class CompType, class SolidThermo>
+inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
+Foam::solidChemistryModel<CompType, SolidThermo>::RR
+(
+    const label i
+) const
+{
+    notImplemented("solidChemistryModel::RR(const label)");
+    return (DimensionedField<scalar, volMesh>::null());
 }
 
 

src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolverType.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -39,23 +39,23 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#define makeSolidChemistrySolverType(SS, Comp, SThermo, GThermo)              \
+#define makeSolidChemistrySolverType(SS, Schem, Comp, SThermo, GThermo)       \
                                                                               \
-    typedef SS<ODESolidChemistryModel<Comp, SThermo, GThermo> >               \
+    typedef SS<Schem<Comp, SThermo, GThermo> >                                \
-        SS##Comp##SThermo##GThermo;                                           \
+        SS##Schem##Comp##SThermo##GThermo;                                    \
                                                                               \
     defineTemplateTypeNameAndDebugWithName                                    \
     (                                                                         \
-        SS##Comp##SThermo##GThermo,                                           \
+        SS##Schem##Comp##SThermo##GThermo,                                    \
-        (#SS"<" + word(Comp::typeName_())                                     \
+        (#SS"<"#Schem"<"#Comp"," + SThermo::typeName() + ","                  \
-      + "," + SThermo::typeName() + "," + GThermo::typeName() + ">").c_str(), \
+      + GThermo::typeName() + ">>").c_str(),                                  \
         0                                                                     \
     );                                                                        \
                                                                               \
     addToRunTimeSelectionTable                                                \
     (                                                                         \
         Comp,                                                                 \
-        SS##Comp##SThermo##GThermo,                                           \
+        SS##Schem##Comp##SThermo##GThermo,                                    \
         fvMesh                                                                \
     );
 

src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolvers.C

@@ -27,8 +27,8 @@ License
 #include "solidThermoPhysicsTypes.H"
 #include "thermoPhysicsTypes.H"
 
-#include "ODESolidChemistryModel.H"
+#include "pyrolysisChemistryModel.H"
-#include "solidChemistryModel.H"
+#include "basicSolidChemistryModel.H"
 
 #include "ode.H"
 
@@ -39,7 +39,8 @@ namespace Foam
     makeSolidChemistrySolverType
     (
         ode,
-        solidChemistryModel,
+        pyrolysisChemistryModel,
+        basicSolidChemistryModel,
         hConstSolidThermoPhysics,
         gasHThermoPhysics
     )
@@ -47,7 +48,8 @@ namespace Foam
     makeSolidChemistrySolverType
     (
         ode,
-        solidChemistryModel,
+        pyrolysisChemistryModel,
+        basicSolidChemistryModel,
         hExponentialSolidThermoPhysics,
         gasHThermoPhysics
     )

src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReaction.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -150,12 +150,18 @@ Foam::string Foam::solidReaction<ReactionThermo>::solidReactionStr
 ) const
 {
     this->reactionStrLeft(reaction);
+    if (glhs().size() > 0)
+    {
         reaction << " + ";
         solidReactionStrLeft(reaction);
+    }
     reaction << " = ";
     this->reactionStrRight(reaction);
+    if (grhs().size() > 0)
+    {
         reaction << " + ";
         solidReactionStrRight(reaction);
+    }
     return reaction.str();
 
 }
@@ -169,8 +175,6 @@ void Foam::solidReaction<ReactionThermo>::solidReactionStrLeft
 {
     for (label i = 0; i < glhs().size(); ++i)
     {
-        reaction << " + ";
-
         if (i > 0)
         {
             reaction << " + ";
@@ -197,8 +201,6 @@ void Foam::solidReaction<ReactionThermo>::solidReactionStrRight
 
     for (label i = 0; i < grhs().size(); ++i)
     {
-        reaction << " + ";
-
         if (i > 0)
         {
             reaction << " + ";