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

applications/test/LduMatrix/LduMatrixTest.C

@@ -0,0 +1,162 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "argList.H"
+#include "Time.H"
+#include "fvMesh.H"
+#include "volFields.H"
+#include "LduMatrix.H"
+#include "diagTensorField.H"
+#include "TPCG.H"
+#include "TPBiCG.H"
+#include "NoPreconditioner.H"
+
+using namespace Foam;
+
+typedef Foam::LduMatrix<vector, diagTensor, scalar>
+    lduVectorMatrix;
+defineNamedTemplateTypeNameAndDebug(lduVectorMatrix, 0);
+
+typedef Foam::DiagonalSolver<vector, diagTensor, scalar>
+    lduVectorDiagonalSolver;
+defineNamedTemplateTypeNameAndDebug(lduVectorDiagonalSolver, 0);
+
+template<>
+const vector lduVectorMatrix::great_(1e15, 1e15, 1e15);
+
+template<>
+const vector lduVectorMatrix::small_(1e-15, 1e-15, 1e-15);
+
+namespace Foam
+{
+    typedef LduMatrix<vector, diagTensor, scalar>::preconditioner
+        lduVectorPreconditioner;
+    defineTemplateRunTimeSelectionTable(lduVectorPreconditioner, symMatrix);
+    defineTemplateRunTimeSelectionTable(lduVectorPreconditioner, asymMatrix);
+
+    typedef LduMatrix<vector, diagTensor, scalar>::smoother
+        lduVectorSmoother;
+    defineTemplateRunTimeSelectionTable(lduVectorSmoother, symMatrix);
+    defineTemplateRunTimeSelectionTable(lduVectorSmoother, asymMatrix);
+
+    typedef LduMatrix<vector, diagTensor, scalar>::solver
+        lduVectorSolver;
+    defineTemplateRunTimeSelectionTable(lduVectorSolver, symMatrix);
+    defineTemplateRunTimeSelectionTable(lduVectorSolver, asymMatrix);
+
+    typedef TPCG<vector, diagTensor, scalar> TPCGVector;
+    defineNamedTemplateTypeNameAndDebug(TPCGVector, 0);
+
+    LduMatrix<vector, diagTensor, scalar>::solver::
+        addsymMatrixConstructorToTable<TPCGVector>
+        addTPCGSymMatrixConstructorToTable_;
+
+    typedef TPBiCG<vector, diagTensor, scalar> TPBiCGVector;
+    defineNamedTemplateTypeNameAndDebug(TPBiCGVector, 0);
+
+    LduMatrix<vector, diagTensor, scalar>::solver::
+        addasymMatrixConstructorToTable<TPBiCGVector>
+        addTPBiCGSymMatrixConstructorToTable_;
+
+    typedef NoPreconditioner<vector, diagTensor, scalar> NoPreconditionerVector;
+    defineNamedTemplateTypeNameAndDebug(NoPreconditionerVector, 0);
+
+    LduMatrix<vector, diagTensor, scalar>::preconditioner::
+        addsymMatrixConstructorToTable<NoPreconditionerVector>
+        addNoPreconditionerSymMatrixConstructorToTable_;
+
+    LduMatrix<vector, diagTensor, scalar>::preconditioner::
+        addasymMatrixConstructorToTable<NoPreconditionerVector>
+        addNoPreconditionerAsymMatrixConstructorToTable_;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Main program:
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+
+    #include "createTime.H"
+    #include "createMesh.H"
+
+    volVectorField psi
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    lduVectorMatrix testMatrix(mesh);
+    testMatrix.diag() = 2*pTraits<diagTensor>::one;
+    testMatrix.source() = pTraits<vector>::one;
+    testMatrix.upper() = 0.1;
+    testMatrix.lower() = -0.1;
+
+    Info<< testMatrix << endl;
+
+    FieldField<Field, scalar> boundaryCoeffs(0);
+    FieldField<Field, scalar> internalCoeffs(0);
+
+    autoPtr<lduVectorMatrix::solver> testMatrixSolver =
+    lduVectorMatrix::solver::New
+    (
+        psi.name(),
+        testMatrix,
+        boundaryCoeffs,
+        internalCoeffs,
+        psi.boundaryField().interfaces(),
+        IStringStream
+        (
+            "PBiCG"
+            "{"
+            "    preconditioner   none;"
+            "    tolerance        (1e-05 1e-05 1e-05);"
+            "    relTol           (0 0 0);"
+            "}"
+        )()
+    );
+
+    lduVectorMatrix::solverPerformance solverPerf =
+        testMatrixSolver->solve(psi);
+
+    solverPerf.print();
+
+    Info<< psi << endl;
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/test/LduMatrix/LduMatrixTest2.C

@@ -0,0 +1,162 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "argList.H"
+#include "Time.H"
+#include "fvMesh.H"
+#include "volFields.H"
+#include "LduMatrix.H"
+#include "tensorField.H"
+#include "TPCG.H"
+#include "TPBiCG.H"
+#include "NoPreconditioner.H"
+
+using namespace Foam;
+
+typedef Foam::LduMatrix<vector, tensor, scalar>
+    lduVectorMatrix;
+defineNamedTemplateTypeNameAndDebug(lduVectorMatrix, 0);
+
+typedef Foam::DiagonalSolver<vector, tensor, scalar>
+    lduVectorDiagonalSolver;
+defineNamedTemplateTypeNameAndDebug(lduVectorDiagonalSolver, 0);
+
+template<>
+const vector lduVectorMatrix::great_(1e15, 1e15, 1e15);
+
+template<>
+const vector lduVectorMatrix::small_(1e-15, 1e-15, 1e-15);
+
+namespace Foam
+{
+    typedef LduMatrix<vector, tensor, scalar>::preconditioner
+        lduVectorPreconditioner;
+    defineTemplateRunTimeSelectionTable(lduVectorPreconditioner, symMatrix);
+    defineTemplateRunTimeSelectionTable(lduVectorPreconditioner, asymMatrix);
+
+    typedef LduMatrix<vector, tensor, scalar>::smoother
+        lduVectorSmoother;
+    defineTemplateRunTimeSelectionTable(lduVectorSmoother, symMatrix);
+    defineTemplateRunTimeSelectionTable(lduVectorSmoother, asymMatrix);
+
+    typedef LduMatrix<vector, tensor, scalar>::solver
+        lduVectorSolver;
+    defineTemplateRunTimeSelectionTable(lduVectorSolver, symMatrix);
+    defineTemplateRunTimeSelectionTable(lduVectorSolver, asymMatrix);
+
+    typedef TPCG<vector, tensor, scalar> TPCGVector;
+    defineNamedTemplateTypeNameAndDebug(TPCGVector, 0);
+
+    LduMatrix<vector, tensor, scalar>::solver::
+        addsymMatrixConstructorToTable<TPCGVector>
+        addTPCGSymMatrixConstructorToTable_;
+
+    typedef TPBiCG<vector, tensor, scalar> TPBiCGVector;
+    defineNamedTemplateTypeNameAndDebug(TPBiCGVector, 0);
+
+    LduMatrix<vector, tensor, scalar>::solver::
+        addasymMatrixConstructorToTable<TPBiCGVector>
+        addTPBiCGSymMatrixConstructorToTable_;
+
+    typedef NoPreconditioner<vector, tensor, scalar> NoPreconditionerVector;
+    defineNamedTemplateTypeNameAndDebug(NoPreconditionerVector, 0);
+
+    LduMatrix<vector, tensor, scalar>::preconditioner::
+        addsymMatrixConstructorToTable<NoPreconditionerVector>
+        addNoPreconditionerSymMatrixConstructorToTable_;
+
+    LduMatrix<vector, tensor, scalar>::preconditioner::
+        addasymMatrixConstructorToTable<NoPreconditionerVector>
+        addNoPreconditionerAsymMatrixConstructorToTable_;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Main program:
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+
+    #include "createTime.H"
+    #include "createMesh.H"
+
+    volVectorField psi
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    lduVectorMatrix testMatrix(mesh);
+    testMatrix.diag() = 2*I;
+    testMatrix.source() = pTraits<vector>::one;
+    testMatrix.upper() = 0.1;
+    testMatrix.lower() = -0.1;
+
+    Info<< testMatrix << endl;
+
+    FieldField<Field, scalar> boundaryCoeffs(0);
+    FieldField<Field, scalar> internalCoeffs(0);
+
+    autoPtr<lduVectorMatrix::solver> testMatrixSolver =
+    lduVectorMatrix::solver::New
+    (
+        psi.name(),
+        testMatrix,
+        //boundaryCoeffs,
+        //internalCoeffs,
+        //psi.boundaryField().interfaces(),
+        IStringStream
+        (
+            "PBiCG"
+            "{"
+            "    preconditioner   none;"
+            "    tolerance        (1e-05 1e-05 1e-05);"
+            "    relTol           (0 0 0);"
+            "}"
+        )()
+    );
+
+    lduVectorMatrix::solverPerformance solverPerf =
+        testMatrixSolver->solve(psi);
+
+    solverPerf.print();
+
+    Info<< psi << endl;
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/test/LduMatrix/LduMatrixTest3.C

@@ -0,0 +1,147 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Application
+    icoFoam
+
+Description
+    Transient solver for incompressible, laminar flow of Newtonian fluids.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "LduMatrix.H"
+#include "diagTensorField.H"
+
+typedef LduMatrix<vector, scalar, scalar> lduVectorMatrix;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        #include "readPISOControls.H"
+        #include "CourantNo.H"
+
+        fvVectorMatrix UEqn
+        (
+            fvm::ddt(U)
+          + fvm::div(phi, U)
+          - fvm::laplacian(nu, U)
+        );
+
+        fvVectorMatrix UEqnp(UEqn == -fvc::grad(p));
+
+        lduVectorMatrix U3Eqnp(mesh);
+        U3Eqnp.diag() = UEqnp.diag();
+        U3Eqnp.upper() = UEqnp.upper();
+        U3Eqnp.lower() = UEqnp.lower();
+        U3Eqnp.source() = UEqnp.source();
+
+        UEqnp.addBoundaryDiag(U3Eqnp.diag(), 0);
+        UEqnp.addBoundarySource(U3Eqnp.source(), false);
+
+        autoPtr<lduVectorMatrix::solver> U3EqnpSolver =
+            lduVectorMatrix::solver::New
+            (
+                U.name(),
+                U3Eqnp,
+                dictionary
+                (
+                    IStringStream
+                    (
+                        "{"
+                        "    solver           PBiCG;"
+                        "    preconditioner   DILU;"
+                        "    tolerance        (1e-13 1e-13 1e-13);"
+                        "    relTol           (0 0 0);"
+                        "}"
+                    )()
+                )
+            );
+
+        U3EqnpSolver->solve(U).print(Info);
+
+        // --- PISO loop
+
+        for (int corr=0; corr<nCorr; corr++)
+        {
+            volScalarField rAU = 1.0/UEqn.A();
+
+            U = rAU*UEqn.H();
+            phi = (fvc::interpolate(U) & mesh.Sf())
+                + fvc::ddtPhiCorr(rAU, U, phi);
+
+            adjustPhi(phi, U, p);
+
+            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+            {
+                fvScalarMatrix pEqn
+                (
+                    fvm::laplacian(rAU, p) == fvc::div(phi)
+                );
+
+                pEqn.setReference(pRefCell, pRefValue);
+                pEqn.solve();
+
+                if (nonOrth == nNonOrthCorr)
+                {
+                    phi -= pEqn.flux();
+                }
+            }
+
+            #include "continuityErrs.H"
+
+            U -= rAU*fvc::grad(p);
+            U.correctBoundaryConditions();
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/test/LduMatrix/Make/files

@@ -0,0 +1,3 @@
+LduMatrixTest3.C
+
+EXE = $(FOAM_USER_APPBIN)/LduMatrixTest

applications/test/LduMatrix/Make/options

@@ -0,0 +1,4 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = -lfiniteVolume

applications/test/LduMatrix/createFields.H

@@ -0,0 +1,55 @@
+    Info<< "Reading transportProperties\n" << endl;
+
+    IOdictionary transportProperties
+    (
+        IOobject
+        (
+            "transportProperties",
+            runTime.constant(),
+            mesh,
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE
+        )
+    );
+
+    dimensionedScalar nu
+    (
+        transportProperties.lookup("nu")
+    );
+
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+
+#   include "createPhi.H"
+
+
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);

src/OpenFOAM/Make/files

@@ -314,6 +314,12 @@ $(algebraicPairGAMGAgglomeration)/algebraicPairGAMGAgglomeration.C
 
 meshes/lduMesh/lduMesh.C
 
+LduMatrix = matrices/LduMatrix
+$(LduMatrix)/LduMatrix/lduMatrices.C
+$(LduMatrix)/Smoothers/lduSmoothers.C
+$(LduMatrix)/Preconditioners/lduPreconditioners.C
+$(LduMatrix)/Solvers/lduSolvers.C
+
 primitiveShapes = meshes/primitiveShapes
 
 $(primitiveShapes)/line/line.C

src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/CyclicLduInterfaceField.C

@@ -0,0 +1,65 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "CyclicLduInterfaceField.H"
+#include "diagTensorField.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    //defineTypeNameAndDebug(CyclicLduInterfaceField, 0);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::CyclicLduInterfaceField<Type>::~CyclicLduInterfaceField()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::CyclicLduInterfaceField<Type>::transformCoupleField
+(
+    Field<Type>& f
+) const
+{
+    if (doTransform())
+    {
+        label sizeby2 = f.size()/2;
+
+        for (label facei=0; facei<sizeby2; facei++)
+        {
+            f[facei] = transform(f[facei], forwardT()[0]);
+            f[facei + sizeby2] = transform(f[facei + sizeby2], forwardT()[0]);
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/CyclicLduInterfaceField.H

@@ -0,0 +1,102 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::CyclicLduInterfaceField
+
+Description
+    Abstract base class for cyclic coupled interfaces.
+
+SourceFiles
+    CyclicLduInterfaceField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef CyclicLduInterfaceField_H
+#define CyclicLduInterfaceField_H
+
+#include "primitiveFieldsFwd.H"
+#include "typeInfo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                  Class CyclicLduInterfaceField Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class CyclicLduInterfaceField
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("CyclicLduInterfaceField");
+
+
+    // Constructors
+
+        //- Construct given coupled patch
+        CyclicLduInterfaceField()
+        {}
+
+
+    // Destructor
+
+        virtual ~CyclicLduInterfaceField();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Is the transform required
+            virtual bool doTransform() const = 0;
+
+            //- Return face transformation tensor
+            virtual const tensorField& forwardT() const = 0;
+
+            //- Return neighbour-cell transformation tensor
+            virtual const tensorField& reverseT() const = 0;
+
+            //- Return rank of component for transform
+            virtual int rank() const = 0;
+
+
+        //- Transform given patch internal field
+        void transformCoupleField(Field<Type>& f) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/LduInterfaceField.C

@@ -0,0 +1,47 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "lduInterfaceField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(lduInterfaceField, 0);
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+lduInterfaceField::~lduInterfaceField()
+{}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/LduInterfaceField.H

@@ -0,0 +1,149 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::LduInterfaceField
+
+Description
+    An abstract base class for implicitly-coupled interface fields
+    e.g. processor and cyclic patch fields.
+
+SourceFiles
+    LduInterfaceField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef LduInterfaceField_H
+#define LduInterfaceField_H
+
+#include "lduInterface.H"
+#include "Field.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class LduInterfaceField Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class LduInterfaceField
+{
+    // Private data
+
+        //- Reference to the coupled patch this field is defined for
+        const lduInterface& interface_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        LduInterfaceField(const LduInterfaceField&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const LduInterfaceField&);
+
+
+public:
+
+        class Amultiplier
+        {
+        public:
+
+            Amultiplier()
+            {}
+
+            virtual ~Amultiplier()
+            {}
+
+            virtual void addAmul
+            (
+                Field<Type>& Apsi,
+                const Field<Type>& psi
+            ) const = 0;
+        };
+
+
+    //- Runtime type information
+    TypeName("LduInterfaceField");
+
+
+    // Constructors
+
+        //- Construct given interface
+        LduInterfaceField(const lduInterface& interface)
+        :
+            interface_(interface)
+        {}
+
+
+    // Destructor
+
+        virtual ~LduInterfaceField();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return the interface
+            const lduInterface& interface() const
+            {
+                return interface_;
+            }
+
+
+        // Coupled interface matrix update
+
+            //- Initialise neighbour matrix update
+            virtual void initInterfaceMatrixUpdate
+            (
+                Field<Type>& Apsi,
+                const Field<Type>& psi,
+                const Amultiplier&,
+                const Pstream::commsTypes commsType
+            ) const
+            {}
+
+            //- Update result field based on interface functionality
+            virtual void updateInterfaceMatrix
+            (
+                Field<Type>& Apsi,
+                const Field<Type>& psi,
+                const Amultiplier&,
+                const Pstream::commsTypes commsType
+            ) const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/LduInterfaceFieldPtrsList.H

@@ -0,0 +1,69 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Type
+    lduInterfaceFieldPtrsList
+
+Description
+    List of coupled interface fields to be used in coupling.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef LduInterfaceFieldPtrsList_H
+#define LduInterfaceFieldPtrsList_H
+
+#include "LduInterfaceField.H"
+#include "UPtrList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class LduInterfaceFieldPtrsList Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class LduInterfaceFieldPtrsList
+:
+    public UPtrList<const LduInterfaceField<Type> >
+{
+public:
+
+    LduInterfaceFieldPtrsList(label size)
+    :
+        UPtrList<const LduInterfaceField<Type> >(size)
+    {}
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/ProcessorLduInterfaceField.C

@@ -0,0 +1,66 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "processorLduInterfaceField.H"
+#include "diagTensorField.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    //defineTypeNameAndDebug(processorLduInterfaceField, 0);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::processorLduInterfaceField<Type>::~processorLduInterfaceField()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::processorLduInterfaceField<Type>::transformCoupleField
+(
+    Field<Type>& f
+) const
+{
+    if (doTransform())
+    {
+        if (forwardT().size() == 1)
+        {
+            transform(f, forwardT()[0], f);
+        }
+        else
+        {
+            transform(f, forwardT(), f);
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/ProcessorLduInterfaceField.H

@@ -0,0 +1,105 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::ProcessorLduInterfaceField
+
+Description
+    Abstract base class for processor coupled interfaces.
+
+SourceFiles
+    ProcessorLduInterfaceField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ProcessorLduInterfaceField_H
+#define ProcessorLduInterfaceField_H
+
+#include "primitiveFieldsFwd.H"
+#include "typeInfo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                  Class ProcessorLduInterfaceField Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class ProcessorLduInterfaceField
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("ProcessorLduInterfaceField");
+
+
+    // Constructors
+
+        //- Construct given coupled patch
+        ProcessorLduInterfaceField()
+        {}
+
+
+    // Destructor
+
+        virtual ~ProcessorLduInterfaceField();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return processor number
+            virtual int myProcNo() const = 0;
+
+            //- Return neigbour processor number
+            virtual int neighbProcNo() const = 0;
+
+            //- Is the transform required
+            virtual bool doTransform() const = 0;
+
+            //- Return face transformation tensor
+            virtual const tensorField& forwardT() const = 0;
+
+            //- Return rank of component for transform
+            virtual int rank() const = 0;
+
+
+        //- Transform given patch component field
+        void transformCoupleField(Field<Type>& f) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.C

@@ -0,0 +1,390 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "lduMatrix.H"
+#include "IOstreams.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+Foam::LduMatrix<Type, DType, LUType>::LduMatrix(const lduMesh& mesh)
+:
+    lduMesh_(mesh),
+    diagPtr_(NULL),
+    upperPtr_(NULL),
+    lowerPtr_(NULL),
+    sourcePtr_(NULL),
+    interfaces_(0),
+    interfacesUpper_(0),
+    interfacesLower_(0)
+{}
+
+
+template<class Type, class DType, class LUType>
+Foam::LduMatrix<Type, DType, LUType>::LduMatrix(const LduMatrix& A)
+:
+    lduMesh_(A.lduMesh_),
+    diagPtr_(NULL),
+    upperPtr_(NULL),
+    lowerPtr_(NULL),
+    sourcePtr_(NULL),
+    interfaces_(0),
+    interfacesUpper_(0),
+    interfacesLower_(0)
+{
+    if (A.diagPtr_)
+    {
+        diagPtr_ = new Field<DType>(*(A.diagPtr_));
+    }
+
+    if (A.upperPtr_)
+    {
+        upperPtr_ = new Field<LUType>(*(A.upperPtr_));
+    }
+
+    if (A.lowerPtr_)
+    {
+        lowerPtr_ = new Field<LUType>(*(A.lowerPtr_));
+    }
+
+    if (A.sourcePtr_)
+    {
+        sourcePtr_ = new Field<Type>(*(A.sourcePtr_));
+    }
+}
+
+
+template<class Type, class DType, class LUType>
+Foam::LduMatrix<Type, DType, LUType>::LduMatrix(LduMatrix& A, bool reUse)
+:
+    lduMesh_(A.lduMesh_),
+    diagPtr_(NULL),
+    upperPtr_(NULL),
+    lowerPtr_(NULL),
+    sourcePtr_(NULL),
+    interfaces_(0),
+    interfacesUpper_(0),
+    interfacesLower_(0)
+{
+    if (reUse)
+    {
+        if (A.diagPtr_)
+        {
+            diagPtr_ = A.diagPtr_;
+            A.diagPtr_ = NULL;
+        }
+
+        if (A.upperPtr_)
+        {
+            upperPtr_ = A.upperPtr_;
+            A.upperPtr_ = NULL;
+        }
+
+        if (A.lowerPtr_)
+        {
+            lowerPtr_ = A.lowerPtr_;
+            A.lowerPtr_ = NULL;
+        }
+
+        if (A.sourcePtr_)
+        {
+            sourcePtr_ = A.sourcePtr_;
+            A.sourcePtr_ = NULL;
+        }
+    }
+    else
+    {
+        if (A.diagPtr_)
+        {
+            diagPtr_ = new Field<DType>(*(A.diagPtr_));
+        }
+
+        if (A.upperPtr_)
+        {
+            upperPtr_ = new Field<LUType>(*(A.upperPtr_));
+        }
+
+        if (A.lowerPtr_)
+        {
+            lowerPtr_ = new Field<LUType>(*(A.lowerPtr_));
+        }
+
+        if (A.sourcePtr_)
+        {
+            sourcePtr_ = new Field<Type>(*(A.sourcePtr_));
+        }
+    }
+}
+
+
+template<class Type, class DType, class LUType>
+Foam::LduMatrix<Type, DType, LUType>::LduMatrix
+(
+    const lduMesh& mesh,
+    Istream& is
+)
+:
+    lduMesh_(mesh),
+    diagPtr_(new Field<DType>(is)),
+    upperPtr_(new Field<LUType>(is)),
+    lowerPtr_(new Field<LUType>(is)),
+    sourcePtr_(new Field<Type>(is)),
+    interfaces_(0),
+    interfacesUpper_(0),
+    interfacesLower_(0)
+{}
+
+
+// * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+Foam::LduMatrix<Type, DType, LUType>::~LduMatrix()
+{
+    if (diagPtr_)
+    {
+        delete diagPtr_;
+    }
+
+    if (upperPtr_)
+    {
+        delete upperPtr_;
+    }
+
+    if (lowerPtr_)
+    {
+        delete lowerPtr_;
+    }
+
+    if (sourcePtr_)
+    {
+        delete sourcePtr_;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+Foam::Field<DType>& Foam::LduMatrix<Type, DType, LUType>::diag()
+{
+    if (!diagPtr_)
+    {
+        diagPtr_ = new Field<DType>(lduAddr().size(), pTraits<DType>::zero);
+    }
+
+    return *diagPtr_;
+}
+
+
+template<class Type, class DType, class LUType>
+Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::upper()
+{
+    if (!upperPtr_)
+    {
+        if (lowerPtr_)
+        {
+            upperPtr_ = new Field<LUType>(*lowerPtr_);
+        }
+        else
+        {
+            upperPtr_ = new Field<LUType>
+            (
+                lduAddr().lowerAddr().size(),
+                pTraits<LUType>::zero
+            );
+        }
+    }
+
+    return *upperPtr_;
+}
+
+
+template<class Type, class DType, class LUType>
+Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::lower()
+{
+    if (!lowerPtr_)
+    {
+        if (upperPtr_)
+        {
+            lowerPtr_ = new Field<LUType>(*upperPtr_);
+        }
+        else
+        {
+            lowerPtr_ = new Field<LUType>
+            (
+                lduAddr().lowerAddr().size(),
+                pTraits<LUType>::zero
+            );
+        }
+    }
+
+    return *lowerPtr_;
+}
+
+
+template<class Type, class DType, class LUType>
+Foam::Field<Type>& Foam::LduMatrix<Type, DType, LUType>::source()
+{
+    if (!sourcePtr_)
+    {
+        sourcePtr_ = new Field<Type>(lduAddr().size(), pTraits<Type>::zero);
+    }
+
+    return *sourcePtr_;
+}
+
+
+template<class Type, class DType, class LUType>
+const Foam::Field<DType>& Foam::LduMatrix<Type, DType, LUType>::diag() const
+{
+    if (!diagPtr_)
+    {
+        FatalErrorIn
+        (
+            "const Field<DType>& LduMatrix<Type, DType, LUType>::diag() const"
+        )   << "diagPtr_ unallocated"
+            << abort(FatalError);
+    }
+
+    return *diagPtr_;
+}
+
+
+template<class Type, class DType, class LUType>
+const Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::upper() const
+{
+    if (!lowerPtr_ && !upperPtr_)
+    {
+        FatalErrorIn
+        (
+            "const Field<LUType>& LduMatrix<Type, DType, LUType>::upper() const"
+        )   << "lowerPtr_ or upperPtr_ unallocated"
+            << abort(FatalError);
+    }
+
+    if (upperPtr_)
+    {
+        return *upperPtr_;
+    }
+    else
+    {
+        return *lowerPtr_;
+    }
+}
+
+
+template<class Type, class DType, class LUType>
+const Foam::Field<LUType>& Foam::LduMatrix<Type, DType, LUType>::lower() const
+{
+    if (!lowerPtr_ && !upperPtr_)
+    {
+        FatalErrorIn
+        (
+            "const Field<LUType>& LduMatrix<Type, DType, LUType>::lower() const"
+        )   << "lowerPtr_ or upperPtr_ unallocated"
+            << abort(FatalError);
+    }
+
+    if (lowerPtr_)
+    {
+        return *lowerPtr_;
+    }
+    else
+    {
+        return *upperPtr_;
+    }
+}
+
+
+template<class Type, class DType, class LUType>
+const Foam::Field<Type>& Foam::LduMatrix<Type, DType, LUType>::source() const
+{
+    if (!sourcePtr_)
+    {
+        FatalErrorIn
+        (
+            "const Field<Type>& LduMatrix<Type, DType, LUType>::source() const"
+        )   << "sourcePtr_ unallocated"
+            << abort(FatalError);
+    }
+
+    return *sourcePtr_;
+}
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const LduMatrix<Type, DType, LUType>& ldum
+)
+{
+    if (ldum.diagPtr_)
+    {
+        os  << "Diagonal = "
+            << *ldum.diagPtr_
+            << endl << endl;
+    }
+
+    if (ldum.upperPtr_)
+    {
+        os  << "Upper triangle = "
+            << *ldum.upperPtr_
+            << endl << endl;
+    }
+
+    if (ldum.lowerPtr_)
+    {
+        os  << "Lower triangle = "
+            << *ldum.lowerPtr_
+            << endl << endl;
+    }
+
+    if (ldum.sourcePtr_)
+    {
+        os  << "Source = "
+            << *ldum.sourcePtr_
+            << endl << endl;
+    }
+
+    os.check("Ostream& operator<<(Ostream&, const LduMatrix&");
+
+    return os;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "LduMatrixOperations.C"
+#include "LduMatrixATmul.C"
+#include "LduMatrixUpdateMatrixInterfaces.C"
+#include "LduMatrixTests.C"
+#include "LduMatrixPreconditioner.C"
+#include "LduMatrixSmoother.C"
+#include "LduMatrixSolver.C"
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H

@@ -0,0 +1,934 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::LduMatrix
+
+Description
+    LduMatrix is a general matrix class in which the coefficients are
+    stored as three arrays, one for the upper triangle, one for the
+    lower triangle and a third for the diagonal.
+
+    Addressing arrays must be supplied for the upper and lower triangles.
+
+Note
+    It might be better if this class were organised as a hierachy starting
+    from an empty matrix, then deriving diagonal, symmetric and asymmetric
+    matrices.
+
+SourceFiles
+    LduMatrixATmul.C
+    LduMatrix.C
+    LduMatrixOperations.C
+    LduMatrixSolver.C
+    LduMatrixPreconditioner.C
+    LduMatrixTests.C
+    LduMatrixUpdateMatrixInterfaces.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef LduMatrix_H
+#define LduMatrix_H
+
+#include "lduMesh.H"
+#include "Field.H"
+#include "FieldField.H"
+#include "LduInterfaceFieldPtrsList.H"
+#include "typeInfo.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of friend functions and operators
+
+template<class Type, class DType, class LUType>
+class LduMatrix;
+
+template<class Type, class DType, class LUType>
+Ostream& operator<<
+(
+    Ostream&,
+    const LduMatrix<Type, DType, LUType>&
+);
+
+
+/*---------------------------------------------------------------------------*\
+                           Class LduMatrix Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type, class DType, class LUType>
+class LduMatrix
+{
+    // private data
+
+        //- LDU mesh reference
+        const lduMesh& lduMesh_;
+
+        //- Diagonal coefficients
+        Field<DType> *diagPtr_;
+
+        //- Off-diagonal coefficients
+        Field<LUType> *upperPtr_, *lowerPtr_;
+
+        //- Source
+        Field<Type> *sourcePtr_;
+
+        //- Field interfaces (processor patches etc.)
+        LduInterfaceFieldPtrsList<Type> interfaces_;
+
+        //- Off-diagonal coefficients for interfaces
+        FieldField<Field, LUType> interfacesUpper_, interfacesLower_;
+
+
+public:
+
+    //- Class returned by the solver
+    //  containing performance statistics
+    class solverPerformance
+    {
+        word   solverName_;
+        word   fieldName_;
+        Type   initialResidual_;
+        Type   finalResidual_;
+        label  noIterations_;
+        bool   converged_;
+        FixedList<bool, pTraits<Type>::nComponents> singular_;
+
+
+    public:
+
+        // Constructors
+
+            solverPerformance()
+            :
+                initialResidual_(pTraits<Type>::zero),
+                finalResidual_(pTraits<Type>::zero),
+                noIterations_(0),
+                converged_(false),
+                singular_(false)
+            {}
+
+
+            solverPerformance
+            (
+                const word&  solverName,
+                const word&  fieldName,
+                const Type&  iRes = pTraits<Type>::zero,
+                const Type&  fRes = pTraits<Type>::zero,
+                const label  nIter = 0,
+                const bool   converged = false,
+                const bool   singular = false
+            )
+            :
+                solverName_(solverName),
+                fieldName_(fieldName),
+                initialResidual_(iRes),
+                finalResidual_(fRes),
+                noIterations_(nIter),
+                converged_(converged),
+                singular_(singular)
+            {}
+
+
+        // Member functions
+
+            //- Return solver name
+            const word& solverName() const
+            {
+                return solverName_;
+            }
+
+            //- Return initial residual
+            const Type& initialResidual() const
+            {
+                return initialResidual_;
+            }
+
+            //- Return initial residual
+            Type& initialResidual()
+            {
+                return initialResidual_;
+            }
+
+
+            //- Return final residual
+            const Type& finalResidual() const
+            {
+                return finalResidual_;
+            }
+
+            //- Return final residual
+            Type& finalResidual()
+            {
+                return finalResidual_;
+            }
+
+
+            //- Return number of iterations
+            label nIterations() const
+            {
+                return noIterations_;
+            }
+
+            //- Return number of iterations
+            label& nIterations()
+            {
+                return noIterations_;
+            }
+
+            //- Check, store and return singularity
+            bool singular(const Type& wApA);
+
+            //- Is the matrix singular?
+            bool singular() const;
+
+            //- Check, store and return convergence
+            bool converged
+            (
+                const Type& tolerance,
+                const Type& relTolerance
+            );
+
+            //- Has the solver converged?
+            bool converged() const
+            {
+                return converged_;
+            }
+
+            //- Print summary of solver performance to the given stream
+            void print(Ostream& os) const;
+    };
+
+
+    //- Abstract base-class for LduMatrix solvers
+    class solver
+    {
+    protected:
+
+        // Protected data
+
+            word fieldName_;
+            const LduMatrix<Type, DType, LUType>& matrix_;
+
+            //- dictionary of controls
+            dictionary controlDict_;
+
+            //- Maximum number of iterations in the solver
+            label maxIter_;
+
+            //- Final convergence tolerance
+            Type tolerance_;
+
+            //- Convergence tolerance relative to the initial
+            Type relTol_;
+
+
+        // Protected Member Functions
+
+            //- Read a control parameter from controlDict
+            template<class T>
+            inline void readControl
+            (
+                const dictionary& controlDict,
+                T& control,
+                const word& controlName
+            );
+
+
+            //- Read the control parameters from the controlDict_
+            virtual void readControls();
+
+
+    public:
+
+        //- Runtime type information
+        virtual const word& type() const = 0;
+
+
+        // Declare run-time constructor selection tables
+
+            declareRunTimeSelectionTable
+            (
+                autoPtr,
+                solver,
+                symMatrix,
+                (
+                    const word& fieldName,
+                    const LduMatrix<Type, DType, LUType>& matrix,
+                    const dictionary& solverDict
+                ),
+                (
+                    fieldName,
+                    matrix,
+                    solverDict
+                )
+            );
+
+            declareRunTimeSelectionTable
+            (
+                autoPtr,
+                solver,
+                asymMatrix,
+                (
+                    const word& fieldName,
+                    const LduMatrix<Type, DType, LUType>& matrix,
+                    const dictionary& solverDict
+                ),
+                (
+                    fieldName,
+                    matrix,
+                    solverDict
+                )
+            );
+
+
+        // Constructors
+
+            solver
+            (
+                const word& fieldName,
+                const LduMatrix<Type, DType, LUType>& matrix,
+                const dictionary& solverDict
+            );
+
+
+        // Selectors
+
+            //- Return a new solver
+            static autoPtr<solver> New
+            (
+                const word& fieldName,
+                const LduMatrix<Type, DType, LUType>& matrix,
+                const dictionary& solverDict
+            );
+
+
+        // Destructor
+
+            virtual ~solver()
+            {}
+
+
+        // Member functions
+
+            // Access
+
+                const word& fieldName() const
+                {
+                    return fieldName_;
+                }
+
+                const LduMatrix<Type, DType, LUType>& matrix() const
+                {
+                    return matrix_;
+                }
+
+
+            //- Read and reset the solver parameters from the given dictionary
+            virtual void read(const dictionary& solverDict);
+
+            virtual solverPerformance solve
+            (
+                Field<Type>& psi
+            ) const = 0;
+
+            //- Return the matrix norm used to normalise the residual for the
+            //  stopping criterion
+            Type normFactor
+            (
+                const Field<Type>& psi,
+                const Field<Type>& Apsi,
+                Field<Type>& tmpField
+            ) const;
+    };
+
+
+    //- Abstract base-class for LduMatrix smoothers
+    class smoother
+    {
+    protected:
+
+        // Protected data
+
+            word fieldName_;
+            const LduMatrix<Type, DType, LUType>& matrix_;
+
+
+    public:
+
+        //- Runtime type information
+        virtual const word& type() const = 0;
+
+
+        // Declare run-time constructor selection tables
+
+            declareRunTimeSelectionTable
+            (
+                autoPtr,
+                smoother,
+                symMatrix,
+                (
+                    const word& fieldName,
+                    const LduMatrix<Type, DType, LUType>& matrix
+                ),
+                (
+                    fieldName,
+                    matrix
+                )
+            );
+
+            declareRunTimeSelectionTable
+            (
+                autoPtr,
+                smoother,
+                asymMatrix,
+                (
+                    const word& fieldName,
+                    const LduMatrix<Type, DType, LUType>& matrix
+                ),
+                (
+                    fieldName,
+                    matrix
+                )
+            );
+
+
+        // Constructors
+
+            smoother
+            (
+                const word& fieldName,
+                const LduMatrix<Type, DType, LUType>& matrix
+            );
+
+
+        // Selectors
+
+            //- Return a new smoother
+            static autoPtr<smoother> New
+            (
+                const word& fieldName,
+                const LduMatrix<Type, DType, LUType>& matrix,
+                const dictionary& smootherDict
+            );
+
+
+        // Destructor
+
+            virtual ~smoother()
+            {}
+
+
+        // Member functions
+
+            // Access
+
+                const word& fieldName() const
+                {
+                    return fieldName_;
+                }
+
+                const LduMatrix<Type, DType, LUType>& matrix() const
+                {
+                    return matrix_;
+                }
+
+
+            //- Smooth the solution for a given number of sweeps
+            virtual void smooth
+            (
+                Field<Type>& psi,
+                const label nSweeps
+            ) const = 0;
+    };
+
+
+    //- Abstract base-class for LduMatrix preconditioners
+    class preconditioner
+    {
+    protected:
+
+        // Protected data
+
+            //- Reference to the base-solver this preconditioner is used with
+            const solver& solver_;
+
+
+    public:
+
+        //- Runtime type information
+        virtual const word& type() const = 0;
+
+
+        // Declare run-time constructor selection tables
+
+            declareRunTimeSelectionTable
+            (
+                autoPtr,
+                preconditioner,
+                symMatrix,
+                (
+                    const solver& sol,
+                    const dictionary& preconditionerDict
+                ),
+                (sol, preconditionerDict)
+            );
+
+            declareRunTimeSelectionTable
+            (
+                autoPtr,
+                preconditioner,
+                asymMatrix,
+                (
+                    const solver& sol,
+                    const dictionary& preconditionerDict
+                ),
+                (sol, preconditionerDict)
+            );
+
+
+        // Constructors
+
+            preconditioner
+            (
+                const solver& sol
+            )
+            :
+                solver_(sol)
+            {}
+
+
+        // Selectors
+
+            //- Return a new preconditioner
+            static autoPtr<preconditioner> New
+            (
+                const solver& sol,
+                const dictionary& preconditionerDict
+            );
+
+
+        // Destructor
+
+            virtual ~preconditioner()
+            {}
+
+
+        // Member functions
+
+            //- Read and reset the preconditioner parameters
+            //  from the given dictionary
+            virtual void read(const dictionary& preconditionerDict)
+            {}
+
+            //- Return wA the preconditioned form of residual rA
+            virtual void precondition
+            (
+                Field<Type>& wA,
+                const Field<Type>& rA
+            ) const = 0;
+
+            //- Return wT the transpose-matrix preconditioned form of
+            //  residual rT.
+            //  This is only required for preconditioning asymmetric matrices.
+            virtual void preconditionT
+            (
+                Field<Type>& wT,
+                const Field<Type>& rT
+            ) const
+            {
+                notImplemented
+                (
+                    type() +"::preconditionT"
+                    "(Field<Type>& wT, const Field<Type>& rT)"
+                );
+            }
+    };
+
+
+    // Static data
+
+        // Declare name of the class and its debug switch
+        ClassName("LduMatrix");
+
+        //- Large Type for the use in solvers
+        static const scalar great_;
+
+        //- Small Type for the use in solvers
+        static const scalar small_;
+
+        //- Very small Type for the use in solvers
+        static const scalar vsmall_;
+
+
+    // Constructors
+
+        //- Construct given an LDU addressed mesh.
+        //  The coefficients are initially empty for subsequent setting.
+        LduMatrix(const lduMesh&);
+
+        //- Construct as copy
+        LduMatrix(const LduMatrix<Type, DType, LUType>&);
+
+        //- Construct as copy or re-use as specified.
+        LduMatrix(LduMatrix<Type, DType, LUType>&, bool reUse);
+
+        //- Construct given an LDU addressed mesh and an Istream
+        //  from which the coefficients are read
+        LduMatrix(const lduMesh&, Istream&);
+
+
+    // Destructor
+
+        ~LduMatrix();
+
+
+    // Member functions
+
+        // Access to addressing
+
+            //- Return the LDU mesh from which the addressing is obtained
+            const lduMesh& mesh() const
+            {
+                return lduMesh_;
+            }
+
+            //- Return the LDU addressing
+            const lduAddressing& lduAddr() const
+            {
+                return lduMesh_.lduAddr();
+            }
+
+            //- Return the patch evaluation schedule
+            const lduSchedule& patchSchedule() const
+            {
+                return lduAddr().patchSchedule();
+            }
+
+            //- Return interfaces
+            const LduInterfaceFieldPtrsList<Type>& interfaces() const
+            {
+                return interfaces_;
+            }
+
+
+        // Access to coefficients
+
+            Field<DType>& diag();
+            Field<LUType>& upper();
+            Field<LUType>& lower();
+            Field<Type>& source();
+
+            FieldField<Field, LUType>& interfacesUpper()
+            {
+                return interfacesUpper_;
+            }
+
+            FieldField<Field, LUType>& interfacesLower()
+            {
+                return interfacesLower_;
+            }
+
+
+            const Field<DType>& diag() const;
+            const Field<LUType>& upper() const;
+            const Field<LUType>& lower() const;
+            const Field<Type>& source() const;
+
+            const FieldField<Field, LUType>& interfacesUpper() const
+            {
+                return interfacesUpper_;
+            }
+
+            const FieldField<Field, LUType>& interfacesLower() const
+            {
+                return interfacesLower_;
+            }
+
+
+            bool hasDiag() const
+            {
+                return (diagPtr_);
+            }
+
+            bool hasUpper() const
+            {
+                return (upperPtr_);
+            }
+
+            bool hasLower() const
+            {
+                return (lowerPtr_);
+            }
+
+            bool hasSource() const
+            {
+                return (sourcePtr_);
+            }
+
+            bool diagonal() const
+            {
+                return (diagPtr_ && !lowerPtr_ && !upperPtr_);
+            }
+
+            bool symmetric() const
+            {
+                return (diagPtr_ && (!lowerPtr_ && upperPtr_));
+            }
+
+            bool asymmetric() const
+            {
+                return (diagPtr_ && lowerPtr_ && upperPtr_);
+            }
+
+
+        // operations
+
+            void sumDiag();
+            void negSumDiag();
+
+            void sumMagOffDiag(Field<LUType>& sumOff) const;
+
+            //- Matrix multiplication
+            void Amul(Field<Type>&, const tmp<Field<Type> >&) const;
+
+            //- Matrix transpose multiplication
+            void Tmul(Field<Type>&, const tmp<Field<Type> >&) const;
+
+
+            //- Sum the coefficients on each row of the matrix
+            void sumA(Field<Type>&) const;
+
+
+            void residual(Field<Type>& rA, const Field<Type>& psi) const;
+
+            tmp<Field<Type> > residual(const Field<Type>& psi) const;
+
+
+            //- Initialise the update of interfaced interfaces
+            //  for matrix operations
+            void initMatrixInterfaces
+            (
+                const FieldField<Field, LUType>& interfaceCoeffs,
+                const Field<Type>& psiif,
+                Field<Type>& result
+            ) const;
+
+            //- Update interfaced interfaces for matrix operations
+            void updateMatrixInterfaces
+            (
+                const FieldField<Field, LUType>& interfaceCoeffs,
+                const Field<Type>& psiif,
+                Field<Type>& result
+            ) const;
+
+
+            tmp<Field<Type> > H(const Field<Type>&) const;
+            tmp<Field<Type> > H(const tmp<Field<Type> >&) const;
+
+            tmp<Field<Type> > faceH(const Field<Type>&) const;
+            tmp<Field<Type> > faceH(const tmp<Field<Type> >&) const;
+
+
+    // Member operators
+
+        void operator=(const LduMatrix<Type, DType, LUType>&);
+
+        void negate();
+
+        void operator+=(const LduMatrix<Type, DType, LUType>&);
+        void operator-=(const LduMatrix<Type, DType, LUType>&);
+
+        void operator*=(const scalarField&);
+        void operator*=(scalar);
+
+
+    // Ostream operator
+
+        friend Ostream& operator<< <Type, DType, LUType>
+        (
+            Ostream&,
+            const LduMatrix<Type, DType, LUType>&
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeLduMatrix(Type, DType, LUType)                                    \
+                                                                              \
+typedef Foam::LduMatrix<Type, DType, LUType>                                  \
+    ldu##Type##DType##LUType##Matrix;                                         \
+                                                                              \
+defineNamedTemplateTypeNameAndDebug(ldu##Type##DType##LUType##Matrix, 0);     \
+                                                                              \
+                                                                              \
+template<>                                                                    \
+const scalar ldu##Type##DType##LUType##Matrix::great_(1e15);                  \
+                                                                              \
+template<>                                                                    \
+const scalar ldu##Type##DType##LUType##Matrix::small_(1e-15);                 \
+                                                                              \
+template<>                                                                    \
+const scalar ldu##Type##DType##LUType##Matrix::vsmall_(VSMALL);               \
+                                                                              \
+                                                                              \
+typedef LduMatrix<Type, DType, LUType>::smoother                              \
+    ldu##Type##DType##LUType##Smoother;                                       \
+                                                                              \
+defineTemplateRunTimeSelectionTable                                           \
+(                                                                             \
+    ldu##Type##DType##LUType##Smoother,                                       \
+    symMatrix                                                                 \
+);                                                                            \
+                                                                              \
+defineTemplateRunTimeSelectionTable                                           \
+(                                                                             \
+    ldu##Type##DType##LUType##Smoother,                                       \
+    asymMatrix                                                                \
+);                                                                            \
+                                                                              \
+                                                                              \
+typedef LduMatrix<Type, DType, LUType>::preconditioner                        \
+    ldu##Type##DType##LUType##Preconditioner;                                 \
+                                                                              \
+defineTemplateRunTimeSelectionTable                                           \
+(                                                                             \
+    ldu##Type##DType##LUType##Preconditioner,                                 \
+    symMatrix                                                                 \
+);                                                                            \
+                                                                              \
+defineTemplateRunTimeSelectionTable                                           \
+(                                                                             \
+    ldu##Type##DType##LUType##Preconditioner,                                 \
+    asymMatrix                                                                \
+);                                                                            \
+                                                                              \
+                                                                              \
+typedef LduMatrix<Type, DType, LUType>::solver                                \
+    ldu##Type##DType##LUType##Solver;                                         \
+                                                                              \
+defineTemplateRunTimeSelectionTable                                           \
+(                                                                             \
+    ldu##Type##DType##LUType##Solver,                                         \
+    symMatrix                                                                 \
+);                                                                            \
+                                                                              \
+defineTemplateRunTimeSelectionTable                                           \
+(                                                                             \
+    ldu##Type##DType##LUType##Solver,                                         \
+    asymMatrix                                                                \
+);
+
+
+#define makeLduPreconditioner(Precon, Type, DType, LUType)                    \
+                                                                              \
+typedef Precon<Type, DType, LUType>                                           \
+    Precon##Type##DType##LUType##Preconditioner;                              \
+defineNamedTemplateTypeNameAndDebug                                           \
+(                                                                             \
+    Precon##Type##DType##LUType##Preconditioner,                              \
+    0                                                                         \
+);
+
+#define makeLduSymPreconditioner(Precon, Type, DType, LUType)                 \
+                                                                              \
+LduMatrix<Type, DType, LUType>::preconditioner::                              \
+addsymMatrixConstructorToTable<Precon##Type##DType##LUType##Preconditioner>   \
+add##Precon##Type##DType##LUType##PreconditionerSymMatrixConstructorToTable_;
+
+#define makeLduAsymPreconditioner(Precon, Type, DType, LUType)                \
+                                                                              \
+LduMatrix<Type, DType, LUType>::preconditioner::                              \
+addasymMatrixConstructorToTable<Precon##Type##DType##LUType##Preconditioner>  \
+add##Precon##Type##DType##LUType##PreconditionerAsymMatrixConstructorToTable_;
+
+
+#define makeLduSmoother(Smoother, Type, DType, LUType)                        \
+                                                                              \
+typedef Smoother<Type, DType, LUType>                                         \
+    Smoother##Type##DType##LUType##Smoother;                                  \
+                                                                              \
+defineNamedTemplateTypeNameAndDebug                                           \
+(                                                                             \
+    Smoother##Type##DType##LUType##Smoother,                                  \
+    0                                                                         \
+);
+
+#define makeLduSymSmoother(Smoother, Type, DType, LUType)                     \
+                                                                              \
+LduMatrix<Type, DType, LUType>::smoother::                                    \
+    addsymMatrixConstructorToTable<Smoother##Type##DType##LUType##Smoother>   \
+    add##Smoother##Type##DType##LUType##SymMatrixConstructorToTable_;
+
+#define makeLduAsymSmoother(Smoother, Type, DType, LUType)                    \
+                                                                              \
+LduMatrix<Type, DType, LUType>::smoother::                                    \
+    addasymMatrixConstructorToTable<Smoother##Type##DType##LUType##Smoother>  \
+    add##Smoother##Type##DType##LUType##AsymMatrixConstructorToTable_;
+
+
+#define makeLduSolver(Solver, Type, DType, LUType)                            \
+                                                                              \
+typedef Solver<Type, DType, LUType>                                           \
+    Solver##Type##DType##LUType##Solver;                                      \
+                                                                              \
+defineNamedTemplateTypeNameAndDebug                                           \
+(                                                                             \
+    Solver##Type##DType##LUType##Solver,                                      \
+    0                                                                         \
+);
+
+#define makeLduSymSolver(Solver, Type, DType, LUType)                         \
+                                                                              \
+LduMatrix<Type, DType, LUType>::solver::                                      \
+    addsymMatrixConstructorToTable<Solver##Type##DType##LUType##Solver>       \
+    add##Solver##Type##DType##LUType##SymMatrixConstructorToTable_;
+
+#define makeLduAsymSolver(Solver, Type, DType, LUType)                        \
+                                                                              \
+LduMatrix<Type, DType, LUType>::solver::                                      \
+    addasymMatrixConstructorToTable<Solver##Type##DType##LUType##Solver>      \
+    add##Solver##Type##DType##LUType##AsymMatrixConstructorToTable_;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "LduMatrixI.H"
+#   include "LduMatrix.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixATmul.C

@@ -0,0 +1,293 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "LduMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<class Type, class LUType>
+class Amultiplier
+:
+    public LduInterfaceField<Type>::Amultiplier
+{
+    const Field<LUType>& A_;
+
+public:
+
+    Amultiplier(const Field<LUType>& A)
+    :
+        A_(A)
+    {}
+
+    virtual ~Amultiplier()
+    {}
+
+    virtual void addAmul(Field<Type>& Apsi, const Field<Type>& psi) const
+    {
+        Apsi += A_*psi;
+    }
+};
+
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+void Foam::LduMatrix<Type, DType, LUType>::Amul
+(
+    Field<Type>& Apsi,
+    const tmp<Field<Type> >& tpsi
+) const
+{
+    Type* __restrict__ ApsiPtr = Apsi.begin();
+
+    const Field<Type>& psi = tpsi();
+    const Type* const __restrict__ psiPtr = psi.begin();
+
+    const DType* const __restrict__ diagPtr = diag().begin();
+
+    const label* const __restrict__ uPtr = lduAddr().upperAddr().begin();
+    const label* const __restrict__ lPtr = lduAddr().lowerAddr().begin();
+
+    const LUType* const __restrict__ upperPtr = upper().begin();
+    const LUType* const __restrict__ lowerPtr = lower().begin();
+
+    // Initialise the update of interfaced interfaces
+    initMatrixInterfaces
+    (
+        interfacesUpper_,
+        psi,
+        Apsi
+    );
+
+    register const label nCells = diag().size();
+    for (register label cell=0; cell<nCells; cell++)
+    {
+        ApsiPtr[cell] = dot(diagPtr[cell], psiPtr[cell]);
+    }
+
+
+    register const label nFaces = upper().size();
+    for (register label face=0; face<nFaces; face++)
+    {
+        ApsiPtr[uPtr[face]] += dot(lowerPtr[face], psiPtr[lPtr[face]]);
+        ApsiPtr[lPtr[face]] += dot(upperPtr[face], psiPtr[uPtr[face]]);
+    }
+
+    // Update interface interfaces
+    updateMatrixInterfaces
+    (
+        interfacesUpper_,
+        psi,
+        Apsi
+    );
+
+    tpsi.clear();
+}
+
+
+template<class Type, class DType, class LUType>
+void Foam::LduMatrix<Type, DType, LUType>::Tmul
+(
+    Field<Type>& Tpsi,
+    const tmp<Field<Type> >& tpsi
+) const
+{
+    Type* __restrict__ TpsiPtr = Tpsi.begin();
+
+    const Field<Type>& psi = tpsi();
+    const Type* const __restrict__ psiPtr = psi.begin();
+
+    const DType* const __restrict__ diagPtr = diag().begin();
+
+    const label* const __restrict__ uPtr = lduAddr().upperAddr().begin();
+    const label* const __restrict__ lPtr = lduAddr().lowerAddr().begin();
+
+    const LUType* const __restrict__ lowerPtr = lower().begin();
+    const LUType* const __restrict__ upperPtr = upper().begin();
+
+    // Initialise the update of interfaced interfaces
+    initMatrixInterfaces
+    (
+        interfacesLower_,
+        psi,
+        Tpsi
+    );
+
+    register const label nCells = diag().size();
+    for (register label cell=0; cell<nCells; cell++)
+    {
+        TpsiPtr[cell] = dot(diagPtr[cell], psiPtr[cell]);
+    }
+
+    register const label nFaces = upper().size();
+    for (register label face=0; face<nFaces; face++)
+    {
+        TpsiPtr[uPtr[face]] += dot(upperPtr[face], psiPtr[lPtr[face]]);
+        TpsiPtr[lPtr[face]] += dot(lowerPtr[face], psiPtr[uPtr[face]]);
+    }
+
+    // Update interface interfaces
+    updateMatrixInterfaces
+    (
+        interfacesLower_,
+        psi,
+        Tpsi
+    );
+
+    tpsi.clear();
+}
+
+
+template<class Type, class DType, class LUType>
+void Foam::LduMatrix<Type, DType, LUType>::sumA
+(
+    Field<Type>& sumA
+) const
+{
+    Type* __restrict__ sumAPtr = sumA.begin();
+
+    const DType* __restrict__ diagPtr = diag().begin();
+
+    const label* __restrict__ uPtr = lduAddr().upperAddr().begin();
+    const label* __restrict__ lPtr = lduAddr().lowerAddr().begin();
+
+    const LUType* __restrict__ lowerPtr = lower().begin();
+    const LUType* __restrict__ upperPtr = upper().begin();
+
+    register const label nCells = diag().size();
+    register const label nFaces = upper().size();
+
+    for (register label cell=0; cell<nCells; cell++)
+    {
+        sumAPtr[cell] = dot(diagPtr[cell], pTraits<Type>::one);
+    }
+
+    for (register label face=0; face<nFaces; face++)
+    {
+        sumAPtr[uPtr[face]] += dot(lowerPtr[face], pTraits<Type>::one);
+        sumAPtr[lPtr[face]] += dot(upperPtr[face], pTraits<Type>::one);
+    }
+
+    // Add the interface internal coefficients to diagonal
+    // and the interface boundary coefficients to the sum-off-diagonal
+    forAll(interfaces_, patchI)
+    {
+        if (interfaces_.set(patchI))
+        {
+            const unallocLabelList& pa = lduAddr().patchAddr(patchI);
+            const Field<LUType>& pCoeffs = interfacesUpper_[patchI];
+
+            forAll(pa, face)
+            {
+                sumAPtr[pa[face]] -= dot(pCoeffs[face], pTraits<Type>::one);
+            }
+        }
+    }
+}
+
+
+template<class Type, class DType, class LUType>
+void Foam::LduMatrix<Type, DType, LUType>::residual
+(
+    Field<Type>& rA,
+    const Field<Type>& psi
+) const
+{
+    Type* __restrict__ rAPtr = rA.begin();
+
+    const Type* const __restrict__ psiPtr = psi.begin();
+    const DType* const __restrict__ diagPtr = diag().begin();
+    const Type* const __restrict__ sourcePtr = source().begin();
+
+    const label* const __restrict__ uPtr = lduAddr().upperAddr().begin();
+    const label* const __restrict__ lPtr = lduAddr().lowerAddr().begin();
+
+    const LUType* const __restrict__ upperPtr = upper().begin();
+    const LUType* const __restrict__ lowerPtr = lower().begin();
+
+    // Parallel boundary initialisation.
+    // Note: there is a change of sign in the coupled
+    // interface update to add the contibution to the r.h.s.
+
+    FieldField<Field, LUType> mBouCoeffs(interfacesUpper_.size());
+
+    forAll(mBouCoeffs, patchi)
+    {
+        if (interfaces_.set(patchi))
+        {
+            mBouCoeffs.set(patchi, -interfacesUpper_[patchi]);
+        }
+    }
+
+    // Initialise the update of interfaced interfaces
+    initMatrixInterfaces
+    (
+        mBouCoeffs,
+        psi,
+        rA
+    );
+
+    register const label nCells = diag().size();
+    for (register label cell=0; cell<nCells; cell++)
+    {
+        rAPtr[cell] = sourcePtr[cell] - dot(diagPtr[cell], psiPtr[cell]);
+    }
+
+
+    register const label nFaces = upper().size();
+    for (register label face=0; face<nFaces; face++)
+    {
+        rAPtr[uPtr[face]] -= dot(lowerPtr[face], psiPtr[lPtr[face]]);
+        rAPtr[lPtr[face]] -= dot(upperPtr[face], psiPtr[uPtr[face]]);
+    }
+
+    // Update interface interfaces
+    updateMatrixInterfaces
+    (
+        mBouCoeffs,
+        psi,
+        rA
+    );
+}
+
+
+template<class Type, class DType, class LUType>
+Foam::tmp<Foam::Field<Type> > Foam::LduMatrix<Type, DType, LUType>::residual
+(
+    const Field<Type>& psi
+) const
+{
+    tmp<Field<Type> > trA(new Field<Type>(psi.size()));
+    residual(trA(), psi);
+    return trA;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixI.H

@@ -0,0 +1,44 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+template<class T>
+inline void Foam::LduMatrix<Type, DType, LUType>::solver::readControl
+(
+    const dictionary& controlDict,
+    T& control,
+    const word& controlName
+)
+{
+    if (controlDict.found(controlName))
+    {
+        controlDict.lookup(controlName) >> control;
+    }
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixOperations.C

@@ -0,0 +1,480 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "lduMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+void Foam::LduMatrix<Type, DType, LUType>::sumDiag()
+{
+    const Field<LUType>& Lower = const_cast<const LduMatrix&>(*this).lower();
+    const Field<LUType>& Upper = const_cast<const LduMatrix&>(*this).upper();
+    Field<DType>& Diag = diag();
+
+    const unallocLabelList& l = lduAddr().lowerAddr();
+    const unallocLabelList& u = lduAddr().upperAddr();
+
+    for (register label face=0; face<l.size(); face++)
+    {
+        Diag[l[face]] += Lower[face];
+        Diag[u[face]] += Upper[face];
+    }
+}
+
+
+template<class Type, class DType, class LUType>
+void Foam::LduMatrix<Type, DType, LUType>::negSumDiag()
+{
+    const Field<LUType>& Lower = const_cast<const LduMatrix&>(*this).lower();
+    const Field<LUType>& Upper = const_cast<const LduMatrix&>(*this).upper();
+    Field<DType>& Diag = diag();
+
+    const unallocLabelList& l = lduAddr().lowerAddr();
+    const unallocLabelList& u = lduAddr().upperAddr();
+
+    for (register label face=0; face<l.size(); face++)
+    {
+        Diag[l[face]] -= Lower[face];
+        Diag[u[face]] -= Upper[face];
+    }
+}
+
+
+template<class Type, class DType, class LUType>
+void Foam::LduMatrix<Type, DType, LUType>::sumMagOffDiag
+(
+    Field<LUType>& sumOff
+) const
+{
+    const Field<LUType>& Lower = const_cast<const LduMatrix&>(*this).lower();
+    const Field<LUType>& Upper = const_cast<const LduMatrix&>(*this).upper();
+
+    const unallocLabelList& l = lduAddr().lowerAddr();
+    const unallocLabelList& u = lduAddr().upperAddr();
+
+    for (register label face = 0; face < l.size(); face++)
+    {
+        sumOff[u[face]] += cmptMag(Lower[face]);
+        sumOff[l[face]] += cmptMag(Upper[face]);
+    }
+}
+
+
+template<class Type, class DType, class LUType>
+Foam::tmp<Foam::Field<Type> >
+Foam::LduMatrix<Type, DType, LUType>::H(const Field<Type>& psi) const
+{
+    tmp<Field<Type> > tHpsi
+    (
+        new Field<Type>(lduAddr().size(), pTraits<Type>::zero)
+    );
+
+    if (lowerPtr_ || upperPtr_)
+    {
+        Field<Type> & Hpsi = tHpsi();
+
+        Type* __restrict__ HpsiPtr = Hpsi.begin();
+
+        const Type* __restrict__ psiPtr = psi.begin();
+
+        const label* __restrict__ uPtr = lduAddr().upperAddr().begin();
+        const label* __restrict__ lPtr = lduAddr().lowerAddr().begin();
+
+        const LUType* __restrict__ lowerPtr = lower().begin();
+        const LUType* __restrict__ upperPtr = upper().begin();
+
+        register const label nFaces = upper().size();
+
+        for (register label face=0; face<nFaces; face++)
+        {
+            HpsiPtr[uPtr[face]] -= lowerPtr[face]*psiPtr[lPtr[face]];
+            HpsiPtr[lPtr[face]] -= upperPtr[face]*psiPtr[uPtr[face]];
+        }
+    }
+
+    return tHpsi;
+}
+
+template<class Type, class DType, class LUType>
+Foam::tmp<Foam::Field<Type> >
+Foam::LduMatrix<Type, DType, LUType>::H(const tmp<Field<Type> >& tpsi) const
+{
+    tmp<Field<Type> > tHpsi(H(tpsi()));
+    tpsi.clear();
+    return tHpsi;
+}
+
+
+template<class Type, class DType, class LUType>
+Foam::tmp<Foam::Field<Type> >
+Foam::LduMatrix<Type, DType, LUType>::faceH(const Field<Type>& psi) const
+{
+    const Field<LUType>& Lower = const_cast<const LduMatrix&>(*this).lower();
+    const Field<LUType>& Upper = const_cast<const LduMatrix&>(*this).upper();
+
+    // Take refereces to addressing
+    const unallocLabelList& l = lduAddr().lowerAddr();
+    const unallocLabelList& u = lduAddr().upperAddr();
+
+    tmp<Field<Type> > tfaceHpsi(new Field<Type> (Lower.size()));
+    Field<Type> & faceHpsi = tfaceHpsi();
+
+    for (register label face=0; face<l.size(); face++)
+    {
+        faceHpsi[face] = Upper[face]*psi[u[face]] - Lower[face]*psi[l[face]];
+    }
+
+    return tfaceHpsi;
+}
+
+
+template<class Type, class DType, class LUType>
+Foam::tmp<Foam::Field<Type> >
+Foam::LduMatrix<Type, DType, LUType>::faceH(const tmp<Field<Type> >& tpsi) const
+{
+    tmp<Field<Type> > tfaceHpsi(faceH(tpsi()));
+    tpsi.clear();
+    return tfaceHpsi;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+void Foam::LduMatrix<Type, DType, LUType>::operator=(const LduMatrix& A)
+{
+    if (this == &A)
+    {
+        FatalErrorIn
+        (
+            "LduMatrix<Type, DType, LUType>::operator=(const LduMatrix&)"
+        )   << "attempted assignment to self"
+            << abort(FatalError);
+    }
+
+    if (A.diagPtr_)
+    {
+        diag() = A.diag();
+    }
+
+    if (A.upperPtr_)
+    {
+        upper() = A.upper();
+    }
+    else if (upperPtr_)
+    {
+        delete upperPtr_;
+        upperPtr_ = NULL;
+    }
+
+    if (A.lowerPtr_)
+    {
+        lower() = A.lower();
+    }
+    else if (lowerPtr_)
+    {
+        delete lowerPtr_;
+        lowerPtr_ = NULL;
+    }
+
+    if (A.sourcePtr_)
+    {
+        source() = A.source();
+    }
+
+    interfacesUpper_ = A.interfacesUpper_;
+    interfacesLower_ = A.interfacesLower_;
+}
+
+
+template<class Type, class DType, class LUType>
+void Foam::LduMatrix<Type, DType, LUType>::negate()
+{
+    if (diagPtr_)
+    {
+        diagPtr_->negate();
+    }
+
+    if (upperPtr_)
+    {
+        upperPtr_->negate();
+    }
+
+    if (lowerPtr_)
+    {
+        lowerPtr_->negate();
+    }
+
+    if (sourcePtr_)
+    {
+        sourcePtr_->negate();
+    }
+
+    negate(interfacesUpper_);
+    negate(interfacesLower_);
+}
+
+
+template<class Type, class DType, class LUType>
+void Foam::LduMatrix<Type, DType, LUType>::operator+=(const LduMatrix& A)
+{
+    if (A.diagPtr_)
+    {
+        diag() += A.diag();
+    }
+
+    if (A.sourcePtr_)
+    {
+        source() += A.source();
+    }
+
+    if (symmetric() && A.symmetric())
+    {
+        upper() += A.upper();
+    }
+    else if (symmetric() && A.asymmetric())
+    {
+        if (upperPtr_)
+        {
+            lower();
+        }
+        else
+        {
+            upper();
+        }
+
+        upper() += A.upper();
+        lower() += A.lower();
+    }
+    else if (asymmetric() && A.symmetric())
+    {
+        if (A.upperPtr_)
+        {
+            lower() += A.upper();
+            upper() += A.upper();
+        }
+        else
+        {
+            lower() += A.lower();
+            upper() += A.lower();
+        }
+
+    }
+    else if (asymmetric() && A.asymmetric())
+    {
+        lower() += A.lower();
+        upper() += A.upper();
+    }
+    else if (diagonal())
+    {
+        if (A.upperPtr_)
+        {
+            upper() = A.upper();
+        }
+
+        if (A.lowerPtr_)
+        {
+            lower() = A.lower();
+        }
+    }
+    else if (A.diagonal())
+    {
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "LduMatrix<Type, DType, LUType>::operator+=(const LduMatrix& A)"
+        )   << "Unknown matrix type combination"
+            << abort(FatalError);
+    }
+
+    interfacesUpper_ += A.interfacesUpper_;
+    interfacesLower_ += A.interfacesLower_;
+}
+
+
+template<class Type, class DType, class LUType>
+void Foam::LduMatrix<Type, DType, LUType>::operator-=(const LduMatrix& A)
+{
+    if (A.diagPtr_)
+    {
+        diag() -= A.diag();
+    }
+
+    if (A.sourcePtr_)
+    {
+        source() -= A.source();
+    }
+
+    if (symmetric() && A.symmetric())
+    {
+        upper() -= A.upper();
+    }
+    else if (symmetric() && A.asymmetric())
+    {
+        if (upperPtr_)
+        {
+            lower();
+        }
+        else
+        {
+            upper();
+        }
+
+        upper() -= A.upper();
+        lower() -= A.lower();
+    }
+    else if (asymmetric() && A.symmetric())
+    {
+        if (A.upperPtr_)
+        {
+            lower() -= A.upper();
+            upper() -= A.upper();
+        }
+        else
+        {
+            lower() -= A.lower();
+            upper() -= A.lower();
+        }
+
+    }
+    else if (asymmetric() && A.asymmetric())
+    {
+        lower() -= A.lower();
+        upper() -= A.upper();
+    }
+    else if (diagonal())
+    {
+        if (A.upperPtr_)
+        {
+            upper() = -A.upper();
+        }
+
+        if (A.lowerPtr_)
+        {
+            lower() = -A.lower();
+        }
+    }
+    else if (A.diagonal())
+    {
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "LduMatrix<Type, DType, LUType>::operator-=(const LduMatrix& A)"
+        )   << "Unknown matrix type combination"
+            << abort(FatalError);
+    }
+
+    interfacesUpper_ -= A.interfacesUpper_;
+    interfacesLower_ -= A.interfacesLower_;
+}
+
+
+template<class Type, class DType, class LUType>
+void Foam::LduMatrix<Type, DType, LUType>::operator*=
+(
+    const scalarField& sf
+)
+{
+    if (diagPtr_)
+    {
+        *diagPtr_ *= sf;
+    }
+
+    if (sourcePtr_)
+    {
+        *sourcePtr_ *= sf;
+    }
+
+    if (upperPtr_)
+    {
+        Field<LUType>& upper = *upperPtr_;
+
+        const unallocLabelList& l = lduAddr().lowerAddr();
+
+        for (register label face=0; face<upper.size(); face++)
+        {
+            upper[face] *= sf[l[face]];
+        }
+    }
+
+    if (lowerPtr_)
+    {
+        Field<LUType>& lower = *lowerPtr_;
+
+        const unallocLabelList& u = lduAddr().upperAddr();
+
+        for (register label face=0; face<lower.size(); face++)
+        {
+            lower[face] *= sf[u[face]];
+        }
+    }
+
+    FatalErrorIn
+    (
+        "LduMatrix<Type, DType, LUType>::operator*=(const scalarField& sf)"
+    )   << "Scaling a matrix by scalarField is not currently supported\n"
+           "because scaling interfacesUpper_ and interfacesLower_ "
+           "require special transfers"
+        << abort(FatalError);
+
+    //interfacesUpper_ *= ;
+    //interfacesLower_ *= sf;
+}
+
+
+template<class Type, class DType, class LUType>
+void Foam::LduMatrix<Type, DType, LUType>::operator*=(scalar s)
+{
+    if (diagPtr_)
+    {
+        *diagPtr_ *= s;
+    }
+
+    if (sourcePtr_)
+    {
+        *sourcePtr_ *= s;
+    }
+
+    if (upperPtr_)
+    {
+        *upperPtr_ *= s;
+    }
+
+    if (lowerPtr_)
+    {
+        *lowerPtr_ *= s;
+    }
+
+    interfacesUpper_ *= s;
+    interfacesLower_ *= s;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixPreconditioner.C

@@ -0,0 +1,115 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "LduMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+Foam::autoPtr<typename Foam::LduMatrix<Type, DType, LUType>::preconditioner>
+Foam::LduMatrix<Type, DType, LUType>::preconditioner::New
+(
+    const solver& sol,
+    const dictionary& preconditionerDict
+)
+{
+    word preconditionerName = preconditionerDict.lookup("preconditioner");
+
+    if (sol.matrix().symmetric())
+    {
+        typename symMatrixConstructorTable::iterator constructorIter =
+            symMatrixConstructorTablePtr_->find(preconditionerName);
+
+        if (constructorIter == symMatrixConstructorTablePtr_->end())
+        {
+            FatalIOErrorIn
+            (
+                "LduMatrix<Type, DType, LUType>::preconditioner::New"
+                "(const solver&, Istream&)",
+                preconditionerDict
+            )   << "Unknown symmetric matrix preconditioner "
+                << preconditionerName << endl << endl
+                << "Valid symmetric matrix preconditioners are :" << endl
+                << symMatrixConstructorTablePtr_->toc()
+                << exit(FatalIOError);
+        }
+
+        return autoPtr<typename LduMatrix<Type, DType, LUType>::preconditioner>
+        (
+            constructorIter()
+            (
+                sol,
+                preconditionerDict
+            )
+        );
+    }
+    else if (sol.matrix().asymmetric())
+    {
+        typename asymMatrixConstructorTable::iterator constructorIter =
+            asymMatrixConstructorTablePtr_->find(preconditionerName);
+
+        if (constructorIter == asymMatrixConstructorTablePtr_->end())
+        {
+            FatalIOErrorIn
+            (
+                "LduMatrix<Type, DType, LUType>::preconditioner::New"
+                "(const solver&, Istream&)",
+                preconditionerDict
+            )   << "Unknown asymmetric matrix preconditioner "
+                << preconditionerName << endl << endl
+                << "Valid asymmetric matrix preconditioners are :" << endl
+                << asymMatrixConstructorTablePtr_->toc()
+                << exit(FatalIOError);
+        }
+
+        return autoPtr<typename LduMatrix<Type, DType, LUType>::preconditioner>
+        (
+            constructorIter()
+            (
+                sol,
+                preconditionerDict
+            )
+        );
+    }
+    else
+    {
+        FatalIOErrorIn
+        (
+            "LduMatrix<Type, DType, LUType>::preconditioner::New"
+            "(const solver&, Istream&)",
+            preconditionerDict
+        )   << "cannot preconditione incomplete matrix, "
+               "no diagonal or off-diagonal coefficient"
+            << exit(FatalIOError);
+
+        return autoPtr<typename LduMatrix<Type, DType, LUType>::preconditioner>
+        (
+            NULL
+        );
+    }
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixSmoother.C

@@ -0,0 +1,120 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "LduMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+Foam::autoPtr<typename Foam::LduMatrix<Type, DType, LUType>::smoother>
+Foam::LduMatrix<Type, DType, LUType>::smoother::New
+(
+    const word& fieldName,
+    const LduMatrix<Type, DType, LUType>& matrix,
+    const dictionary& smootherDict
+)
+{
+    word smootherName = smootherDict.lookup("smoother");
+
+    if (matrix.symmetric())
+    {
+        typename symMatrixConstructorTable::iterator constructorIter =
+            symMatrixConstructorTablePtr_->find(smootherName);
+
+        if (constructorIter == symMatrixConstructorTablePtr_->end())
+        {
+            FatalIOErrorIn
+            (
+                "LduMatrix<Type, DType, LUType>::smoother::New", smootherDict
+            )   << "Unknown symmetric matrix smoother " << smootherName
+                << endl << endl
+                << "Valid symmetric matrix smoothers are :" << endl
+                << symMatrixConstructorTablePtr_->toc()
+                << exit(FatalIOError);
+        }
+
+        return autoPtr<typename LduMatrix<Type, DType, LUType>::smoother>
+        (
+            constructorIter()
+            (
+                fieldName,
+                matrix
+            )
+        );
+    }
+    else if (matrix.asymmetric())
+    {
+        typename asymMatrixConstructorTable::iterator constructorIter =
+            asymMatrixConstructorTablePtr_->find(smootherName);
+
+        if (constructorIter == asymMatrixConstructorTablePtr_->end())
+        {
+            FatalIOErrorIn
+            (
+                "LduMatrix<Type, DType, LUType>::smoother::New", smootherDict
+            )   << "Unknown asymmetric matrix smoother " << smootherName
+                << endl << endl
+                << "Valid asymmetric matrix smoothers are :" << endl
+                << asymMatrixConstructorTablePtr_->toc()
+                << exit(FatalIOError);
+        }
+
+        return autoPtr<typename LduMatrix<Type, DType, LUType>::smoother>
+        (
+            constructorIter()
+            (
+                fieldName,
+                matrix
+            )
+        );
+    }
+    else
+    {
+        FatalIOErrorIn
+        (
+            "LduMatrix<Type, DType, LUType>::smoother::New", smootherDict
+        )   << "cannot solve incomplete matrix, no off-diagonal coefficients"
+            << exit(FatalIOError);
+
+        return autoPtr<typename LduMatrix<Type, DType, LUType>::smoother>(NULL);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+Foam::LduMatrix<Type, DType, LUType>::smoother::smoother
+(
+    const word& fieldName,
+    const LduMatrix<Type, DType, LUType>& matrix
+)
+:
+    fieldName_(fieldName),
+    matrix_(matrix)
+{}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixSolver.C

@@ -0,0 +1,190 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "LduMatrix.H"
+#include "DiagonalSolver.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+Foam::autoPtr<typename Foam::LduMatrix<Type, DType, LUType>::solver>
+Foam::LduMatrix<Type, DType, LUType>::solver::New
+(
+    const word& fieldName,
+    const LduMatrix<Type, DType, LUType>& matrix,
+    const dictionary& solverDict
+)
+{
+    word solverName = solverDict.lookup("solver");
+
+    if (matrix.diagonal())
+    {
+        return autoPtr<typename LduMatrix<Type, DType, LUType>::solver>
+        (
+            new DiagonalSolver<Type, DType, LUType>
+            (
+                fieldName,
+                matrix,
+                solverDict
+            )
+        );
+    }
+    else if (matrix.symmetric())
+    {
+        typename symMatrixConstructorTable::iterator constructorIter =
+            symMatrixConstructorTablePtr_->find(solverName);
+
+        if (constructorIter == symMatrixConstructorTablePtr_->end())
+        {
+            FatalIOErrorIn
+            (
+                "LduMatrix<Type, DType, LUType>::solver::New", solverDict
+            )   << "Unknown symmetric matrix solver " << solverName
+                << endl << endl
+                << "Valid symmetric matrix solvers are :" << endl
+                << symMatrixConstructorTablePtr_->toc()
+                << exit(FatalIOError);
+        }
+
+        return autoPtr<typename LduMatrix<Type, DType, LUType>::solver>
+        (
+            constructorIter()
+            (
+                fieldName,
+                matrix,
+                solverDict
+            )
+        );
+    }
+    else if (matrix.asymmetric())
+    {
+        typename asymMatrixConstructorTable::iterator constructorIter =
+            asymMatrixConstructorTablePtr_->find(solverName);
+
+        if (constructorIter == asymMatrixConstructorTablePtr_->end())
+        {
+            FatalIOErrorIn
+            (
+                "LduMatrix<Type, DType, LUType>::solver::New", solverDict
+            )   << "Unknown asymmetric matrix solver " << solverName
+                << endl << endl
+                << "Valid asymmetric matrix solvers are :" << endl
+                << asymMatrixConstructorTablePtr_->toc()
+                << exit(FatalIOError);
+        }
+
+        return autoPtr<typename LduMatrix<Type, DType, LUType>::solver>
+        (
+            constructorIter()
+            (
+                fieldName,
+                matrix,
+                solverDict
+            )
+        );
+    }
+    else
+    {
+        FatalIOErrorIn
+        (
+            "LduMatrix<Type, DType, LUType>::solver::New", solverDict
+        )   << "cannot solve incomplete matrix, "
+               "no diagonal or off-diagonal coefficient"
+            << exit(FatalIOError);
+
+        return autoPtr<typename LduMatrix<Type, DType, LUType>::solver>(NULL);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+Foam::LduMatrix<Type, DType, LUType>::solver::solver
+(
+    const word& fieldName,
+    const LduMatrix<Type, DType, LUType>& matrix,
+    const dictionary& solverDict
+)
+:
+    fieldName_(fieldName),
+    matrix_(matrix),
+
+    controlDict_(solverDict),
+
+    maxIter_(1000),
+    tolerance_(1e-6*pTraits<Type>::one),
+    relTol_(pTraits<Type>::zero)
+{
+    readControls();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+void Foam::LduMatrix<Type, DType, LUType>::solver::readControls()
+{
+    readControl(controlDict_, maxIter_, "maxIter");
+    readControl(controlDict_, tolerance_, "tolerance");
+    readControl(controlDict_, relTol_, "relTol");
+}
+
+
+template<class Type, class DType, class LUType>
+void Foam::LduMatrix<Type, DType, LUType>::solver::read
+(
+    const dictionary& solverDict
+)
+{
+    controlDict_ = solverDict;
+    readControls();
+}
+
+
+template<class Type, class DType, class LUType>
+Type Foam::LduMatrix<Type, DType, LUType>::solver::normFactor
+(
+    const Field<Type>& psi,
+    const Field<Type>& Apsi,
+    Field<Type>& tmpField
+) const
+{
+    // --- Calculate A dot reference value of psi
+    matrix_.sumA(tmpField);
+    cmptMultiply(tmpField, tmpField, gAverage(psi));
+
+    return stabilise
+    (
+        gSum(cmptMag(Apsi - tmpField) + cmptMag(matrix_.source() - tmpField)),
+        matrix_.small_
+    );
+
+    // At convergence this simpler method is equivalent to the above
+    // return stabilise(2*gSumCmptMag(matrix_.source()), matrix_.small_);
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixTests.C

@@ -0,0 +1,118 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "LduMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+bool Foam::LduMatrix<Type, DType, LUType>::solverPerformance::singular
+(
+    const Type& wApA
+)
+{
+    for(direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
+    {
+        singular_[cmpt] = component(wApA, cmpt) < vsmall_;
+    }
+
+    return singular();
+}
+
+
+template<class Type, class DType, class LUType>
+bool Foam::LduMatrix<Type, DType, LUType>::solverPerformance::singular() const
+{
+    for(direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
+    {
+        if (!singular_[cmpt]) return false;
+    }
+
+    return true;
+}
+
+
+template<class Type, class DType, class LUType>
+bool Foam::LduMatrix<Type, DType, LUType>::solverPerformance::converged
+(
+    const Type& Tolerance,
+    const Type& RelTolerance
+)
+{
+    if (debug >= 2)
+    {
+        Info<< solverName_
+            << ":  Iteration " << noIterations_
+            << " residual = " << finalResidual_
+            << endl;
+    }
+
+    if
+    (
+        finalResidual_ < Tolerance
+     || (
+            RelTolerance > small_*pTraits<Type>::one
+         && finalResidual_ < cmptMultiply(RelTolerance, initialResidual_)
+        )
+    )
+    {
+        converged_ = true;
+    }
+    else
+    {
+        converged_ = false;
+    }
+
+    return converged_;
+}
+
+
+template<class Type, class DType, class LUType>
+void Foam::LduMatrix<Type, DType, LUType>::solverPerformance::print
+(
+    Ostream& os
+) const
+{
+    for(direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
+    {
+        os  << solverName_ << ":  Solving for "
+            << word(fieldName_ + pTraits<Type>::componentNames[cmpt]);
+
+        if (singular_[cmpt])
+        {
+            os  << ":  solution singularity" << endl;
+        }
+        else
+        {
+            os  << ", Initial residual = " << component(initialResidual_, cmpt)
+                << ", Final residual = " << component(finalResidual_, cmpt)
+                << ", No Iterations " << noIterations_
+                << endl;
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixUpdateMatrixInterfaces.C

@@ -0,0 +1,194 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "LduMatrix.H"
+#include "LduInterfaceField.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+void Foam::LduMatrix<Type, DType, LUType>::initMatrixInterfaces
+(
+    const FieldField<Field, LUType>& interfaceCoeffs,
+    const Field<Type>& psiif,
+    Field<Type>& result
+) const
+{
+    if
+    (
+        Pstream::defaultCommsType == Pstream::blocking
+     || Pstream::defaultCommsType == Pstream::nonBlocking
+    )
+    {
+        forAll (interfaces_, interfaceI)
+        {
+            if (interfaces_.set(interfaceI))
+            {
+                interfaces_[interfaceI].initInterfaceMatrixUpdate
+                (
+                    result,
+                    psiif,
+                    Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
+                    Pstream::defaultCommsType
+                );
+            }
+        }
+    }
+    else if (Pstream::defaultCommsType == Pstream::scheduled)
+    {
+        const lduSchedule& patchSchedule = this->patchSchedule();
+
+        // Loop over the "global" patches are on the list of interfaces but
+        // beyond the end of the schedule which only handles "normal" patches
+        for
+        (
+            label interfaceI=patchSchedule.size()/2;
+            interfaceI<interfaces_.size();
+            interfaceI++
+        )
+        {
+            if (interfaces_.set(interfaceI))
+            {
+                interfaces_[interfaceI].initInterfaceMatrixUpdate
+                (
+                    result,
+                    psiif,
+                    Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
+                    Pstream::blocking
+                );
+            }
+        }
+    }
+    else
+    {
+        FatalErrorIn("LduMatrix<Type, DType, LUType>::initMatrixInterfaces")
+            << "Unsuported communications type "
+            << Pstream::commsTypeNames[Pstream::defaultCommsType]
+            << exit(FatalError);
+    }
+}
+
+
+template<class Type, class DType, class LUType>
+void Foam::LduMatrix<Type, DType, LUType>::updateMatrixInterfaces
+(
+    const FieldField<Field, LUType>& interfaceCoeffs,
+    const Field<Type>& psiif,
+    Field<Type>& result
+) const
+{
+    if
+    (
+        Pstream::defaultCommsType == Pstream::blocking
+     || Pstream::defaultCommsType == Pstream::nonBlocking
+    )
+    {
+        // Block until all sends/receives have been finished
+        if (Pstream::defaultCommsType == Pstream::nonBlocking)
+        {
+            IPstream::waitRequests();
+            OPstream::waitRequests();
+        }
+
+        forAll (interfaces_, interfaceI)
+        {
+            if (interfaces_.set(interfaceI))
+            {
+                interfaces_[interfaceI].updateInterfaceMatrix
+                (
+                    result,
+                    psiif,
+                    Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
+                    Pstream::defaultCommsType
+                );
+            }
+        }
+    }
+    else if (Pstream::defaultCommsType == Pstream::scheduled)
+    {
+        const lduSchedule& patchSchedule = this->patchSchedule();
+
+        // Loop over all the "normal" interfaces relating to standard patches
+        forAll (patchSchedule, i)
+        {
+            label interfaceI = patchSchedule[i].patch;
+
+            if (interfaces_.set(interfaceI))
+            {
+                if (patchSchedule[i].init)
+                {
+                    interfaces_[interfaceI].initInterfaceMatrixUpdate
+                    (
+                        result,
+                        psiif,
+                        Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
+                        Pstream::scheduled
+                    );
+                }
+                else
+                {
+                    interfaces_[interfaceI].updateInterfaceMatrix
+                    (
+                        result,
+                        psiif,
+                        Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
+                        Pstream::scheduled
+                    );
+                }
+            }
+        }
+
+        // Loop over the "global" patches are on the list of interfaces but
+        // beyond the end of the schedule which only handles "normal" patches
+        for
+        (
+            label interfaceI=patchSchedule.size()/2;
+            interfaceI<interfaces_.size();
+            interfaceI++
+        )
+        {
+            if (interfaces_.set(interfaceI))
+            {
+                interfaces_[interfaceI].updateInterfaceMatrix
+                (
+                    result,
+                    psiif,
+                    Amultiplier<Type, LUType>(interfaceCoeffs[interfaceI]),
+                    Pstream::blocking
+                );
+            }
+        }
+    }
+    else
+    {
+        FatalErrorIn("LduMatrix<Type, DType, LUType>::updateMatrixInterfaces")
+            << "Unsuported communications type "
+            << Pstream::commsTypeNames[Pstream::defaultCommsType]
+            << exit(FatalError);
+    }
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/LduMatrix/lduMatrices.C

@@ -0,0 +1,39 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "LduMatrix.H"
+#include "fieldTypes.H"
+
+namespace Foam
+{
+    makeLduMatrix(scalar, scalar, scalar);
+    makeLduMatrix(vector, scalar, scalar);
+    makeLduMatrix(sphericalTensor, scalar, scalar);
+    makeLduMatrix(symmTensor, scalar, scalar);
+    makeLduMatrix(tensor, scalar, scalar);
+};
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/Preconditioners/DILUPreconditioner/TDILUPreconditioner.C

@@ -0,0 +1,179 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "TDILUPreconditioner.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+Foam::TDILUPreconditioner<Type, DType, LUType>::TDILUPreconditioner
+(
+    const typename LduMatrix<Type, DType, LUType>::solver& sol,
+    const dictionary&
+)
+:
+    LduMatrix<Type, DType, LUType>::preconditioner(sol),
+    rD_(sol.matrix().diag())
+{
+    calcInvD(rD_, sol.matrix());
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+void Foam::TDILUPreconditioner<Type, DType, LUType>::calcInvD
+(
+    Field<DType>& rD,
+    const LduMatrix<Type, DType, LUType>& matrix
+)
+{
+    DType* __restrict__ rDPtr = rD.begin();
+
+    const label* const __restrict__ uPtr = matrix.lduAddr().upperAddr().begin();
+    const label* const __restrict__ lPtr = matrix.lduAddr().lowerAddr().begin();
+
+    const LUType* const __restrict__ upperPtr = matrix.upper().begin();
+    const LUType* const __restrict__ lowerPtr = matrix.lower().begin();
+
+    register label nFaces = matrix.upper().size();
+    for (register label face=0; face<nFaces; face++)
+    {
+        rDPtr[uPtr[face]] -=
+            dot(dot(upperPtr[face], lowerPtr[face]), inv(rDPtr[lPtr[face]]));
+    }
+
+
+    // Calculate the reciprocal of the preconditioned diagonal
+    register label nCells = rD.size();
+
+    for (register label cell=0; cell<nCells; cell++)
+    {
+        rDPtr[cell] = inv(rDPtr[cell]);
+    }
+}
+
+
+template<class Type, class DType, class LUType>
+void Foam::TDILUPreconditioner<Type, DType, LUType>::precondition
+(
+    Field<Type>& wA,
+    const Field<Type>& rA
+) const
+{
+    Type* __restrict__ wAPtr = wA.begin();
+    const Type* __restrict__ rAPtr = rA.begin();
+    const DType* __restrict__ rDPtr = rD_.begin();
+
+    const label* const __restrict__ uPtr =
+        this->solver_.matrix().lduAddr().upperAddr().begin();
+    const label* const __restrict__ lPtr =
+        this->solver_.matrix().lduAddr().lowerAddr().begin();
+    const label* const __restrict__ losortPtr =
+        this->solver_.matrix().lduAddr().losortAddr().begin();
+
+    const LUType* const __restrict__ upperPtr =
+        this->solver_.matrix().upper().begin();
+    const LUType* const __restrict__ lowerPtr =
+        this->solver_.matrix().lower().begin();
+
+    register label nCells = wA.size();
+    register label nFaces = this->solver_.matrix().upper().size();
+    register label nFacesM1 = nFaces - 1;
+
+    for (register label cell=0; cell<nCells; cell++)
+    {
+        wAPtr[cell] = dot(rDPtr[cell], rAPtr[cell]);
+    }
+
+
+    register label sface;
+
+    for (register label face=0; face<nFaces; face++)
+    {
+        sface = losortPtr[face];
+        wAPtr[uPtr[sface]] -=
+            dot(rDPtr[uPtr[sface]], dot(lowerPtr[sface], wAPtr[lPtr[sface]]));
+    }
+
+    for (register label face=nFacesM1; face>=0; face--)
+    {
+        wAPtr[lPtr[face]] -=
+            dot(rDPtr[lPtr[face]], dot(upperPtr[face], wAPtr[uPtr[face]]));
+    }
+}
+
+
+template<class Type, class DType, class LUType>
+void Foam::TDILUPreconditioner<Type, DType, LUType>::preconditionT
+(
+    Field<Type>& wT,
+    const Field<Type>& rT
+) const
+{
+    Type* __restrict__ wTPtr = wT.begin();
+    const Type* __restrict__ rTPtr = rT.begin();
+    const DType* __restrict__ rDPtr = rD_.begin();
+
+    const label* const __restrict__ uPtr =
+        this->solver_.matrix().lduAddr().upperAddr().begin();
+    const label* const __restrict__ lPtr =
+        this->solver_.matrix().lduAddr().lowerAddr().begin();
+    const label* const __restrict__ losortPtr =
+        this->solver_.matrix().lduAddr().losortAddr().begin();
+
+    const LUType* const __restrict__ upperPtr =
+        this->solver_.matrix().upper().begin();
+    const LUType* const __restrict__ lowerPtr =
+        this->solver_.matrix().lower().begin();
+
+    register label nCells = wT.size();
+    register label nFaces = this->solver_.matrix().upper().size();
+    register label nFacesM1 = nFaces - 1;
+
+    for (register label cell=0; cell<nCells; cell++)
+    {
+        wTPtr[cell] = dot(rDPtr[cell], rTPtr[cell]);
+    }
+
+    for (register label face=0; face<nFaces; face++)
+    {
+        wTPtr[uPtr[face]] -=
+            dot(rDPtr[uPtr[face]], dot(upperPtr[face], wTPtr[lPtr[face]]));
+    }
+
+
+    register label sface;
+
+    for (register label face=nFacesM1; face>=0; face--)
+    {
+        sface = losortPtr[face];
+        wTPtr[lPtr[sface]] -=
+            dot(rDPtr[lPtr[sface]], dot(lowerPtr[sface], wTPtr[uPtr[sface]]));
+    }
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/Preconditioners/DILUPreconditioner/TDILUPreconditioner.H

@@ -0,0 +1,126 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::TDILUPreconditioner
+
+Description
+    Simplified diagonal-based incomplete LU preconditioner for asymmetric
+    matrices.
+
+    The inverse (reciprocal for scalar) of the preconditioned diagonal is
+    calculated and stored.
+
+SourceFiles
+    TDILUPreconditioner.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef TDILUPreconditioner_H
+#define TDILUPreconditioner_H
+
+#include "LduMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class TDILUPreconditioner Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type, class DType, class LUType>
+class TDILUPreconditioner
+:
+    public LduMatrix<Type, DType, LUType>::preconditioner
+{
+    // Private data
+
+        //- The inverse (reciprocal for scalar) preconditioned diagonal
+        Field<DType> rD_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("DILU");
+
+
+    // Constructors
+
+        //- Construct from matrix components and preconditioner data dictionary
+        TDILUPreconditioner
+        (
+            const typename LduMatrix<Type, DType, LUType>::solver& sol,
+            const dictionary& preconditionerDict
+        );
+
+
+    // Destructor
+
+        virtual ~TDILUPreconditioner()
+        {}
+
+
+    // Member Functions
+
+        //- Calculate the reciprocal of the preconditioned diagonal
+        static void calcInvD
+        (
+            Field<DType>& rD,
+            const LduMatrix<Type, DType, LUType>& matrix
+        );
+
+        //- Return wA the preconditioned form of residual rA
+        virtual void precondition
+        (
+            Field<Type>& wA,
+            const Field<Type>& rA
+        ) const;
+
+        //- Return wT the transpose-matrix preconditioned form of
+        //  residual rT.
+        virtual void preconditionT
+        (
+            Field<Type>& wT,
+            const Field<Type>& rT
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "TDILUPreconditioner.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/Preconditioners/DiagonalPreconditioner/DiagonalPreconditioner.C

@@ -0,0 +1,80 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "DiagonalPreconditioner.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+Foam::DiagonalPreconditioner<Type, DType, LUType>::DiagonalPreconditioner
+(
+    const typename LduMatrix<Type, DType, LUType>::solver& sol,
+    const dictionary&
+)
+:
+    LduMatrix<Type, DType, LUType>::preconditioner(sol),
+    rD(sol.matrix().diag().size())
+{
+    DType* __restrict__ rDPtr = rD.begin();
+    const DType* __restrict__ DPtr = this->solver_.matrix().diag().begin();
+
+    register label nCells = rD.size();
+
+    // Generate inverse (reciprocal for scalar) diagonal
+    for (register label cell=0; cell<nCells; cell++)
+    {
+        rDPtr[cell] = inv(DPtr[cell]);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+void Foam::DiagonalPreconditioner<Type, DType, LUType>::read(const dictionary&)
+{}
+
+
+template<class Type, class DType, class LUType>
+void Foam::DiagonalPreconditioner<Type, DType, LUType>::precondition
+(
+    Field<Type>& wA,
+    const Field<Type>& rA
+) const
+{
+    Type* __restrict__ wAPtr = wA.begin();
+    const Type* __restrict__ rAPtr = rA.begin();
+    const DType* __restrict__ rDPtr = rD.begin();
+
+    register label nCells = wA.size();
+
+    for (register label cell=0; cell<nCells; cell++)
+    {
+        wAPtr[cell] = dot(rDPtr[cell], rAPtr[cell]);
+    }
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/Preconditioners/DiagonalPreconditioner/DiagonalPreconditioner.H

@@ -0,0 +1,134 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::DiagonalPreconditioner
+
+Description
+    Diagonal preconditioner for both symmetric and asymmetric matrices.
+
+    The inverse (reciprocal for scalar) of the diagonal is calculated and
+    stored.
+
+SourceFiles
+    DiagonalPreconditioner.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef DiagonalPreconditioner_H
+#define DiagonalPreconditioner_H
+
+#include "LduMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class DiagonalPreconditioner Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type, class DType, class LUType>
+class DiagonalPreconditioner
+:
+    public LduMatrix<Type, DType, LUType>::preconditioner
+{
+    // Private data
+
+        //- The inverse (reciprocal for scalar) diagonal
+        Field<DType> rD;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        DiagonalPreconditioner(const DiagonalPreconditioner&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const DiagonalPreconditioner&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("diagonal");
+
+
+    // Constructors
+
+        //- Construct from matrix components and preconditioner data dictionary
+        DiagonalPreconditioner
+        (
+            const typename LduMatrix<Type, DType, LUType>::solver& sol,
+            const dictionary& preconditionerDict
+        );
+
+
+    // Destructor
+
+        virtual ~DiagonalPreconditioner()
+        {}
+
+
+    // Member Functions
+
+        //- Read and reset the preconditioner parameters from the given
+        //  dictionary
+        virtual void read(const dictionary& preconditionerDict);
+
+        //- Return wA the preconditioned form of residual rA
+        virtual void precondition
+        (
+            Field<Type>& wA,
+            const Field<Type>& rA
+        ) const;
+
+        //- Return wT the transpose-matrix preconditioned form of
+        //  residual rT.
+        virtual void preconditionT
+        (
+            Field<Type>& wT,
+            const Field<Type>& rT
+        ) const
+        {
+            return(precondition(wT, rT));
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "DiagonalPreconditioner.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/Preconditioners/NoPreconditioner/NoPreconditioner.C

@@ -0,0 +1,59 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "NoPreconditioner.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+Foam::NoPreconditioner<Type, DType, LUType>::NoPreconditioner
+(
+    const typename LduMatrix<Type, DType, LUType>::solver& sol,
+    const dictionary&
+)
+:
+    LduMatrix<Type, DType, LUType>::preconditioner(sol)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+void Foam::NoPreconditioner<Type, DType, LUType>::read(const dictionary&)
+{}
+
+
+template<class Type, class DType, class LUType>
+void Foam::NoPreconditioner<Type, DType, LUType>::precondition
+(
+    Field<Type>& wA,
+    const Field<Type>& rA
+) const
+{
+    wA = rA;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/Preconditioners/NoPreconditioner/NoPreconditioner.H

@@ -0,0 +1,125 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::NoPreconditioner
+
+Description
+    Null preconditioner for both symmetric and asymmetric matrices.
+
+SourceFiles
+    NoPreconditioner.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef NoPreconditioner_H
+#define NoPreconditioner_H
+
+#include "LduMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class NoPreconditioner Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type, class DType, class LUType>
+class NoPreconditioner
+:
+    public LduMatrix<Type, DType, LUType>::preconditioner
+{
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        NoPreconditioner(const NoPreconditioner&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const NoPreconditioner&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("none");
+
+
+    // Constructors
+
+        //- Construct from matrix components and preconditioner data dictionary
+        NoPreconditioner
+        (
+            const typename LduMatrix<Type, DType, LUType>::solver& sol,
+            const dictionary& preconditionerDict
+        );
+
+
+    // Destructor
+
+        virtual ~NoPreconditioner()
+        {}
+
+
+    // Member Functions
+
+        //- Read and reset the preconditioner parameters from the given
+        //  dictionary
+        virtual void read(const dictionary& preconditionerDict);
+
+        //- Return wA the preconditioned form of residual rA
+        virtual void precondition
+        (
+            Field<Type>& wA,
+            const Field<Type>& rA
+        ) const;
+
+        //- Return wT the transpose-matrix preconditioned form of
+        //  residual rT.
+        virtual void preconditionT
+        (
+            Field<Type>& wT,
+            const Field<Type>& rT
+        ) const
+        {
+            return(precondition(wT, rT));
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "NoPreconditioner.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/Preconditioners/lduPreconditioners.C

@@ -0,0 +1,54 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "NoPreconditioner.H"
+#include "DiagonalPreconditioner.H"
+#include "TDILUPreconditioner.H"
+#include "fieldTypes.H"
+
+#define makeLduPreconditioners(Type, DType, LUType)                           \
+                                                                              \
+    makeLduPreconditioner(NoPreconditioner, Type, DType, LUType);             \
+    makeLduSymPreconditioner(NoPreconditioner, Type, DType, LUType);          \
+    makeLduAsymPreconditioner(NoPreconditioner, Type, DType, LUType);         \
+                                                                              \
+    makeLduPreconditioner(DiagonalPreconditioner, Type, DType, LUType);       \
+    makeLduSymPreconditioner(DiagonalPreconditioner, Type, DType, LUType);    \
+    makeLduAsymPreconditioner(DiagonalPreconditioner, Type, DType, LUType);   \
+                                                                              \
+    makeLduPreconditioner(TDILUPreconditioner, Type, DType, LUType);          \
+    makeLduAsymPreconditioner(TDILUPreconditioner, Type, DType, LUType);
+
+namespace Foam
+{
+    makeLduPreconditioners(scalar, scalar, scalar);
+    makeLduPreconditioners(vector, scalar, scalar);
+    makeLduPreconditioners(sphericalTensor, scalar, scalar);
+    makeLduPreconditioners(symmTensor, scalar, scalar);
+    makeLduPreconditioners(tensor, scalar, scalar);
+};
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/Smoothers/GaussSeidel/TGaussSeidelSmoother.C

@@ -0,0 +1,167 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "TGaussSeidelSmoother.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+Foam::TGaussSeidelSmoother<Type, DType, LUType>::TGaussSeidelSmoother
+(
+    const word& fieldName,
+    const LduMatrix<Type, DType, LUType>& matrix
+)
+:
+    LduMatrix<Type, DType, LUType>::smoother
+    (
+        fieldName,
+        matrix
+    ),
+    rD_(matrix.diag().size())
+{
+    register const label nCells = matrix.diag().size();
+    register const DType* const __restrict__ diagPtr = matrix.diag().begin();
+    register DType* __restrict__ rDPtr = rD_.begin();
+
+    for (register label cellI=0; cellI<nCells; cellI++)
+    {
+        rDPtr[cellI] = inv(diagPtr[cellI]);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+void Foam::TGaussSeidelSmoother<Type, DType, LUType>::smooth
+(
+    const word& fieldName_,
+    Field<Type>& psi,
+    const LduMatrix<Type, DType, LUType>& matrix_,
+    const Field<DType>& rD_,
+    const label nSweeps
+)
+{
+    register Type* __restrict__ psiPtr = psi.begin();
+
+    register const label nCells = psi.size();
+
+    Field<Type> bPrime(nCells);
+    register Type* __restrict__ bPrimePtr = bPrime.begin();
+
+    register const DType* const __restrict__ rDPtr = rD_.begin();
+
+    register const LUType* const __restrict__ upperPtr =
+        matrix_.upper().begin();
+
+    register const LUType* const __restrict__ lowerPtr =
+        matrix_.lower().begin();
+
+    register const label* const __restrict__ uPtr =
+        matrix_.lduAddr().upperAddr().begin();
+
+    register const label* const __restrict__ ownStartPtr =
+        matrix_.lduAddr().ownerStartAddr().begin();
+
+
+    // Parallel boundary initialisation.  The parallel boundary is treated
+    // as an effective jacobi interface in the boundary.
+    // Note: there is a change of sign in the coupled
+    // interface update to add the contibution to the r.h.s.
+
+    FieldField<Field, LUType> mBouCoeffs(matrix_.interfacesUpper().size());
+
+    forAll(mBouCoeffs, patchi)
+    {
+        if (matrix_.interfaces().set(patchi))
+        {
+            mBouCoeffs.set(patchi, -matrix_.interfacesUpper()[patchi]);
+        }
+    }
+
+    for (label sweep=0; sweep<nSweeps; sweep++)
+    {
+        bPrime = matrix_.source();
+
+        matrix_.initMatrixInterfaces
+        (
+            mBouCoeffs,
+            psi,
+            bPrime
+        );
+
+        matrix_.updateMatrixInterfaces
+        (
+            mBouCoeffs,
+            psi,
+            bPrime
+        );
+
+        Type curPsi;
+        register label fStart;
+        register label fEnd = ownStartPtr[0];
+
+        for (register label cellI=0; cellI<nCells; cellI++)
+        {
+            // Start and end of this row
+            fStart = fEnd;
+            fEnd = ownStartPtr[cellI + 1];
+
+            // Get the accumulated neighbour side
+            curPsi = bPrimePtr[cellI];
+
+            // Accumulate the owner product side
+            for (register label curFace=fStart; curFace<fEnd; curFace++)
+            {
+                curPsi -= dot(upperPtr[curFace], psiPtr[uPtr[curFace]]);
+            }
+
+            // Finish current psi
+            curPsi = dot(rDPtr[cellI], curPsi);
+
+            // Distribute the neighbour side using current psi
+            for (register label curFace=fStart; curFace<fEnd; curFace++)
+            {
+                bPrimePtr[uPtr[curFace]] -= dot(lowerPtr[curFace], curPsi);
+            }
+
+            psiPtr[cellI] = curPsi;
+        }
+    }
+}
+
+
+template<class Type, class DType, class LUType>
+void Foam::TGaussSeidelSmoother<Type, DType, LUType>::smooth
+(
+    Field<Type>& psi,
+    const label nSweeps
+) const
+{
+    smooth(this->fieldName_, psi, this->matrix_, rD_, nSweeps);
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/Smoothers/GaussSeidel/TGaussSeidelSmoother.H

@@ -0,0 +1,112 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::TGaussSeidelSmoother
+
+Description
+    Foam::TGaussSeidelSmoother
+
+SourceFiles
+    TGaussSeidelSmoother.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef TGaussSeidelSmoother_H
+#define TGaussSeidelSmoother_H
+
+#include "LduMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class TGaussSeidelSmoother Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type, class DType, class LUType>
+class TGaussSeidelSmoother
+:
+    public LduMatrix<Type, DType, LUType>::smoother
+{
+    // Private data
+
+        //- The inverse (reciprocal for scalars) diagonal
+        Field<DType> rD_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("GaussSeidel");
+
+
+    // Constructors
+
+        //- Construct from components
+        TGaussSeidelSmoother
+        (
+            const word& fieldName,
+            const LduMatrix<Type, DType, LUType>& matrix
+        );
+
+
+    // Member Functions
+
+        //- Smooth for the given number of sweeps
+        static void smooth
+        (
+            const word& fieldName,
+            Field<Type>& psi,
+            const LduMatrix<Type, DType, LUType>& matrix,
+            const Field<DType>& rD,
+            const label nSweeps
+        );
+
+
+        //- Smooth the solution for a given number of sweeps
+        virtual void smooth
+        (
+            Field<Type>& psi,
+            const label nSweeps
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "TGaussSeidelSmoother.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/Smoothers/lduSmoothers.C

@@ -0,0 +1,45 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "TGaussSeidelSmoother.H"
+#include "fieldTypes.H"
+
+#define makeLduSmoothers(Type, DType, LUType)                                 \
+                                                                              \
+    makeLduSmoother(TGaussSeidelSmoother, Type, DType, LUType);               \
+    makeLduSymSmoother(TGaussSeidelSmoother, Type, DType, LUType);            \
+    makeLduAsymSmoother(TGaussSeidelSmoother, Type, DType, LUType);
+
+namespace Foam
+{
+    makeLduSmoothers(scalar, scalar, scalar);
+    makeLduSmoothers(vector, scalar, scalar);
+    makeLduSmoothers(sphericalTensor, scalar, scalar);
+    makeLduSmoothers(symmTensor, scalar, scalar);
+    makeLduSmoothers(tensor, scalar, scalar);
+};
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.C

@@ -0,0 +1,79 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "DiagonalSolver.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+Foam::DiagonalSolver<Type, DType, LUType>::DiagonalSolver
+(
+    const word& fieldName,
+    const LduMatrix<Type, DType, LUType>& matrix,
+    const dictionary& solverDict
+)
+:
+    LduMatrix<Type, DType, LUType>::solver
+    (
+        fieldName,
+        matrix,
+        solverDict
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+void Foam::DiagonalSolver<Type, DType, LUType>::read
+(
+    const dictionary&
+)
+{}
+
+
+template<class Type, class DType, class LUType>
+typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance
+Foam::DiagonalSolver<Type, DType, LUType>::solve
+(
+    Field<Type>& psi
+) const
+{
+    psi = this->matrix_.source()/this->matrix_.diag();
+
+    return typename LduMatrix<Type, DType, LUType>::solverPerformance
+    (
+        typeName,
+        this->fieldName_,
+        pTraits<Type>::zero,
+        pTraits<Type>::zero,
+        0,
+        true,
+        false
+    );
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.H

@@ -0,0 +1,107 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::DiagonalSolver
+
+Description
+    Foam::DiagonalSolver
+
+SourceFiles
+    DiagonalSolver.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef DiagonalSolver_H
+#define DiagonalSolver_H
+
+#include "LduMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class DiagonalSolver Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type, class DType, class LUType>
+class DiagonalSolver
+:
+    public LduMatrix<Type, DType, LUType>::solver
+{
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        DiagonalSolver(const DiagonalSolver&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const DiagonalSolver&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("diagonal");
+
+
+    // Constructors
+
+        //- Construct from matrix
+        DiagonalSolver
+        (
+            const word& fieldName,
+            const LduMatrix<Type, DType, LUType>& matrix,
+            const dictionary& solverDict
+        );
+
+
+    // Member Functions
+
+        //- Read and reset the solver parameters from the given dictionary
+        void read(const dictionary& solverDict);
+
+        //- Solve the matrix with this solver
+        typename LduMatrix<Type, DType, LUType>::solverPerformance solve
+        (
+            Field<Type>& psi
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "DiagonalSolver.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCG/TPBiCG.C

@@ -0,0 +1,194 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "TPBiCG.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+Foam::TPBiCG<Type, DType, LUType>::TPBiCG
+(
+    const word& fieldName,
+    const LduMatrix<Type, DType, LUType>& matrix,
+    const dictionary& solverDict
+)
+:
+    LduMatrix<Type, DType, LUType>::solver
+    (
+        fieldName,
+        matrix,
+        solverDict
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance
+Foam::TPBiCG<Type, DType, LUType>::solve(Field<Type>& psi) const
+{
+    word preconditionerName(this->controlDict_.lookup("preconditioner"));
+
+    // --- Setup class containing solver performance data
+    typename LduMatrix<Type, DType, LUType>::solverPerformance solverPerf
+    (
+        preconditionerName + typeName,
+        this->fieldName_
+    );
+
+    register label nCells = psi.size();
+
+    Type* __restrict__ psiPtr = psi.begin();
+
+    Field<Type> pA(nCells);
+    Type* __restrict__ pAPtr = pA.begin();
+
+    Field<Type> pT(nCells, pTraits<Type>::zero);
+    Type* __restrict__ pTPtr = pT.begin();
+
+    Field<Type> wA(nCells);
+    Type* __restrict__ wAPtr = wA.begin();
+
+    Field<Type> wT(nCells);
+    Type* __restrict__ wTPtr = wT.begin();
+
+    Type wArT = this->matrix_.great_*pTraits<Type>::one;
+    Type wArTold = wArT;
+
+    // --- Calculate A.psi and T.psi
+    this->matrix_.Amul(wA, psi);
+    this->matrix_.Tmul(wT, psi);
+
+    // --- Calculate initial residual and transpose residual fields
+    Field<Type> rA(this->matrix_.source() - wA);
+    Field<Type> rT(this->matrix_.source() - wT);
+    Type* __restrict__ rAPtr = rA.begin();
+    Type* __restrict__ rTPtr = rT.begin();
+
+    // --- Calculate normalisation factor
+    Type normFactor = this->normFactor(psi, wA, pA);
+
+    if (LduMatrix<Type, DType, LUType>::debug >= 2)
+    {
+        Info<< "   Normalisation factor = " << normFactor << endl;
+    }
+
+    // --- Calculate normalised residual norm
+    solverPerf.initialResidual() = cmptDivide(gSumCmptMag(rA), normFactor);
+    solverPerf.finalResidual() = solverPerf.initialResidual();
+
+    // --- Check convergence, solve if not converged
+    if (!solverPerf.converged(this->tolerance_, this->relTol_))
+    {
+        // --- Select and construct the preconditioner
+        autoPtr<typename LduMatrix<Type, DType, LUType>::preconditioner>
+        preconPtr = LduMatrix<Type, DType, LUType>::preconditioner::New
+        (
+            *this,
+            this->controlDict_
+        );
+
+        // --- Solver iteration
+        do
+        {
+            // --- Store previous wArT
+            wArTold = wArT;
+
+            // --- Precondition residuals
+            preconPtr->precondition(wA, rA);
+            preconPtr->preconditionT(wT, rT);
+
+            // --- Update search directions:
+            wArT = gSumCmptProd(wA, rT);
+
+            if (solverPerf.nIterations() == 0)
+            {
+                for (register label cell=0; cell<nCells; cell++)
+                {
+                    pAPtr[cell] = wAPtr[cell];
+                    pTPtr[cell] = wTPtr[cell];
+                }
+            }
+            else
+            {
+                Type beta = cmptDivide
+                (
+                    wArT,
+                    stabilise(wArTold, this->matrix_.vsmall_)
+                );
+
+                for (register label cell=0; cell<nCells; cell++)
+                {
+                    pAPtr[cell] = wAPtr[cell] + cmptMultiply(beta, pAPtr[cell]);
+                    pTPtr[cell] = wTPtr[cell] + cmptMultiply(beta, pTPtr[cell]);
+                }
+            }
+
+
+            // --- Update preconditioned residuals
+            this->matrix_.Amul(wA, pA);
+
+            this->matrix_.Tmul(wT, pT);
+
+            Type wApT = gSumCmptProd(wA, pT);
+
+            // --- Test for singularity
+            if (solverPerf.singular(cmptDivide(cmptMag(wApT), normFactor)))
+            {
+                break;
+            }
+
+
+            // --- Update solution and residual:
+
+            Type alpha = cmptDivide
+            (
+                wArT,
+                stabilise(wApT, this->matrix_.vsmall_)
+            );
+
+            for (register label cell=0; cell<nCells; cell++)
+            {
+                psiPtr[cell] += cmptMultiply(alpha, pAPtr[cell]);
+                rAPtr[cell] -= cmptMultiply(alpha, wAPtr[cell]);
+                rTPtr[cell] -= cmptMultiply(alpha, wTPtr[cell]);
+            }
+
+            solverPerf.finalResidual() =
+                cmptDivide(gSumCmptMag(rA), normFactor);
+
+        } while
+        (
+            solverPerf.nIterations()++ < this->maxIter_
+        && !(solverPerf.converged(this->tolerance_, this->relTol_))
+        );
+    }
+
+    return solverPerf;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCG/TPBiCG.H

@@ -0,0 +1,111 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::TPBiCG
+
+Description
+    Preconditioned bi-conjugate gradient solver for asymmetric lduMatrices
+    using a run-time selectable preconditiioner.
+
+SourceFiles
+    TPBiCG.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef TPBiCG_H
+#define TPBiCG_H
+
+#include "LduMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class TPBiCG Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type, class DType, class LUType>
+class TPBiCG
+:
+    public LduMatrix<Type, DType, LUType>::solver
+{
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        TPBiCG(const TPBiCG&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const TPBiCG&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("PBiCG");
+
+
+    // Constructors
+
+        //- Construct from matrix components and solver data dictionary
+        TPBiCG
+        (
+            const word& fieldName,
+            const LduMatrix<Type, DType, LUType>& matrix,
+            const dictionary& solverDict
+        );
+
+
+    // Destructor
+
+        virtual ~TPBiCG()
+        {}
+
+
+    // Member Functions
+
+        //- Solve the matrix with this solver
+        virtual typename LduMatrix<Type, DType, LUType>::solverPerformance solve
+        (
+            Field<Type>& psi
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "TPBiCG.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCGScalarAlpha/PBiCGScalarAlpha.C

@@ -0,0 +1,190 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "TPBiCG.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+Foam::TPBiCG<Type, DType, LUType>::TPBiCG
+(
+    const word& fieldName,
+    const LduMatrix<Type, DType, LUType>& matrix,
+    const dictionary& solverDict
+)
+:
+    LduMatrix<Type, DType, LUType>::solver
+    (
+        fieldName,
+        matrix,
+        solverDict
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance
+Foam::TPBiCG<Type, DType, LUType>::solve
+(
+    Field<Type>& psi
+) const
+{
+    word preconditionerName(this->controlDict_.lookup("preconditioner"));
+
+    // --- Setup class containing solver performance data
+    typename LduMatrix<Type, DType, LUType>::solverPerformance solverPerf
+    (
+        preconditionerName + typeName,
+        this->fieldName_
+    );
+
+    register label nCells = psi.size();
+
+    Type* __restrict__ psiPtr = psi.begin();
+
+    Field<Type> pA(nCells);
+    Type* __restrict__ pAPtr = pA.begin();
+
+    Field<Type> pT(nCells, pTraits<Type>::zero);
+    Type* __restrict__ pTPtr = pT.begin();
+
+    Field<Type> wA(nCells);
+    Type* __restrict__ wAPtr = wA.begin();
+
+    Field<Type> wT(nCells);
+    Type* __restrict__ wTPtr = wT.begin();
+
+    scalar wArT = 1e15; //this->matrix_.great_;
+    scalar wArTold = wArT;
+
+    // --- Calculate A.psi and T.psi
+    this->matrix_.Amul(wA, psi);
+    this->matrix_.Tmul(wT, psi);
+
+    // --- Calculate initial residual and transpose residual fields
+    Field<Type> rA(this->matrix_.source() - wA);
+    Field<Type> rT(this->matrix_.source() - wT);
+    Type* __restrict__ rAPtr = rA.begin();
+    Type* __restrict__ rTPtr = rT.begin();
+
+    // --- Calculate normalisation factor
+    Type normFactor = this->normFactor(psi, wA, pA);
+
+    if (LduMatrix<Type, DType, LUType>::debug >= 2)
+    {
+        Info<< "   Normalisation factor = " << normFactor << endl;
+    }
+
+    // --- Calculate normalised residual norm
+    solverPerf.initialResidual() = cmptDivide(gSumCmptMag(rA), normFactor);
+    solverPerf.finalResidual() = solverPerf.initialResidual();
+
+    // --- Check convergence, solve if not converged
+    if (!solverPerf.converged(this->tolerance_, this->relTol_))
+    {
+        // --- Select and construct the preconditioner
+        autoPtr<typename LduMatrix<Type, DType, LUType>::preconditioner>
+        preconPtr = LduMatrix<Type, DType, LUType>::preconditioner::New
+        (
+            *this,
+            this->controlDict_
+        );
+
+        // --- Solver iteration
+        do
+        {
+            // --- Store previous wArT
+            wArTold = wArT;
+
+            // --- Precondition residuals
+            preconPtr->precondition(wA, rA);
+            preconPtr->preconditionT(wT, rT);
+
+            // --- Update search directions:
+            //wArT = gSumProd(wA, rT);
+            wArT = sum(wA & rT);
+
+            if (solverPerf.nIterations() == 0)
+            {
+                for (register label cell=0; cell<nCells; cell++)
+                {
+                    pAPtr[cell] = wAPtr[cell];
+                    pTPtr[cell] = wTPtr[cell];
+                }
+            }
+            else
+            {
+                scalar beta = cmptDivide(wArT, wArTold);
+
+                for (register label cell=0; cell<nCells; cell++)
+                {
+                    pAPtr[cell] = wAPtr[cell] + (beta* pAPtr[cell]);
+                    pTPtr[cell] = wTPtr[cell] + (beta* pTPtr[cell]);
+                }
+            }
+
+
+            // --- Update preconditioned residuals
+            this->matrix_.Amul(wA, pA);
+            this->matrix_.Tmul(wT, pT);
+
+            scalar wApT = sum(wA & pT);
+
+
+            // --- Test for singularity
+            //if
+            //(
+            //    solverPerf.checkSingularity(ratio(mag(wApT)normFactor))
+            //) break;
+
+
+            // --- Update solution and residual:
+
+            scalar alpha = cmptDivide(wArT, wApT);
+
+            for (register label cell=0; cell<nCells; cell++)
+            {
+                psiPtr[cell] += (alpha* pAPtr[cell]);
+                rAPtr[cell] -= (alpha* wAPtr[cell]);
+                rTPtr[cell] -= (alpha* wTPtr[cell]);
+            }
+
+            solverPerf.finalResidual() =
+                cmptDivide(gSumCmptMag(rA), normFactor);
+
+        } while
+        (
+            solverPerf.nIterations()++ < this->maxIter_
+        && !(solverPerf.converged(this->tolerance_, this->relTol_))
+        );
+    }
+
+    return solverPerf;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/Solvers/PBiCGScalarAlpha/PBiCGScalarAlpha.H

@@ -0,0 +1,111 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::TPBiCG
+
+Description
+    Preconditioned bi-conjugate gradient solver for asymmetric lduMatrices
+    using a run-time selectable preconditiioner.
+
+SourceFiles
+    TPBiCG.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef TPBiCG_H
+#define TPBiCG_H
+
+#include "LduMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class TPBiCG Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type, class DType, class LUType>
+class TPBiCG
+:
+    public LduMatrix<Type, DType, LUType>::solver
+{
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        TPBiCG(const TPBiCG&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const TPBiCG&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("PBiCG");
+
+
+    // Constructors
+
+        //- Construct from matrix components and solver data dictionary
+        TPBiCG
+        (
+            const word& fieldName,
+            const LduMatrix<Type, DType, LUType>& matrix,
+            const dictionary& solverDict
+        );
+
+
+    // Destructor
+
+        virtual ~TPBiCG()
+        {}
+
+
+    // Member Functions
+
+        //- Solve the matrix with this solver
+        virtual typename LduMatrix<Type, DType, LUType>::solverPerformance solve
+        (
+            Field<Type>& psi
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "PBiCGScalarAlpha.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/Solvers/PCG/TPCG.C

@@ -0,0 +1,180 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "TPCG.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+Foam::TPCG<Type, DType, LUType>::TPCG
+(
+    const word& fieldName,
+    const LduMatrix<Type, DType, LUType>& matrix,
+    const dictionary& solverDict
+)
+:
+    LduMatrix<Type, DType, LUType>::solver
+    (
+        fieldName,
+        matrix,
+        solverDict
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance
+Foam::TPCG<Type, DType, LUType>::solve(Field<Type>& psi) const
+{
+    word preconditionerName(this->controlDict_.lookup("preconditioner"));
+
+    // --- Setup class containing solver performance data
+    typename LduMatrix<Type, DType, LUType>::solverPerformance solverPerf
+    (
+        preconditionerName + typeName,
+        this->fieldName_
+    );
+
+    register label nCells = psi.size();
+
+    Type* __restrict__ psiPtr = psi.begin();
+
+    Field<Type> pA(nCells);
+    Type* __restrict__ pAPtr = pA.begin();
+
+    Field<Type> wA(nCells);
+    Type* __restrict__ wAPtr = wA.begin();
+
+    Type wArA = this->matrix_.great_*pTraits<Type>::one;
+    Type wArAold = wArA;
+
+    // --- Calculate A.psi
+    this->matrix_.Amul(wA, psi);
+
+    // --- Calculate initial residual field
+    Field<Type> rA(this->matrix_.source() - wA);
+    Type* __restrict__ rAPtr = rA.begin();
+
+    // --- Calculate normalisation factor
+    Type normFactor = this->normFactor(psi, wA, pA);
+
+    if (LduMatrix<Type, DType, LUType>::debug >= 2)
+    {
+        Info<< "   Normalisation factor = " << normFactor << endl;
+    }
+
+    // --- Calculate normalised residual norm
+    solverPerf.initialResidual() = cmptDivide(gSumCmptMag(rA), normFactor);
+    solverPerf.finalResidual() = solverPerf.initialResidual();
+
+    // --- Check convergence, solve if not converged
+    if (!solverPerf.converged(this->tolerance_, this->relTol_))
+    {
+        // --- Select and construct the preconditioner
+        autoPtr<typename LduMatrix<Type, DType, LUType>::preconditioner>
+        preconPtr = LduMatrix<Type, DType, LUType>::preconditioner::New
+        (
+            *this,
+            this->controlDict_
+        );
+
+        // --- Solver iteration
+        do
+        {
+            // --- Store previous wArA
+            wArAold = wArA;
+
+            // --- Precondition residual
+            preconPtr->precondition(wA, rA);
+
+            // --- Update search directions:
+            wArA = gSumCmptProd(wA, rA);
+
+            if (solverPerf.nIterations() == 0)
+            {
+                for (register label cell=0; cell<nCells; cell++)
+                {
+                    pAPtr[cell] = wAPtr[cell];
+                }
+            }
+            else
+            {
+                Type beta = cmptDivide
+                (
+                    wArA,
+                    stabilise(wArAold, this->matrix_.vsmall_)
+                );
+
+                for (register label cell=0; cell<nCells; cell++)
+                {
+                    pAPtr[cell] = wAPtr[cell] + cmptMultiply(beta, pAPtr[cell]);
+                }
+            }
+
+
+            // --- Update preconditioned residual
+            this->matrix_.Amul(wA, pA);
+
+            Type wApA = gSumCmptProd(wA, pA);
+
+
+            // --- Test for singularity
+            if (solverPerf.singular(cmptDivide(cmptMag(wApA), normFactor)))
+            {
+                break;
+            }
+
+
+            // --- Update solution and residual:
+
+            Type alpha = cmptDivide
+            (
+                wArA,
+                stabilise(wApA, this->matrix_.vsmall_)
+            );
+
+            for (register label cell=0; cell<nCells; cell++)
+            {
+                psiPtr[cell] += cmptMultiply(alpha, pAPtr[cell]);
+                rAPtr[cell] -= cmptMultiply(alpha, wAPtr[cell]);
+            }
+
+            solverPerf.finalResidual() =
+                cmptDivide(gSumCmptMag(rA), normFactor);
+
+        } while
+        (
+            solverPerf.nIterations()++ < this->maxIter_
+        && !(solverPerf.converged(this->tolerance_, this->relTol_))
+        );
+    }
+
+    return solverPerf;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/Solvers/PCG/TPCG.H

@@ -0,0 +1,111 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::TPCG
+
+Description
+    Preconditioned conjugate gradient solver for symmetric lduMatrices
+    using a run-time selectable preconditiioner.
+
+SourceFiles
+    TPCG.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef TPCG_H
+#define TPCG_H
+
+#include "LduMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class TPCG Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type, class DType, class LUType>
+class TPCG
+:
+    public LduMatrix<Type, DType, LUType>::solver
+{
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        TPCG(const TPCG&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const TPCG&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("PCG");
+
+
+    // Constructors
+
+        //- Construct from matrix components and solver data dictionary
+        TPCG
+        (
+            const word& fieldName,
+            const LduMatrix<Type, DType, LUType>& matrix,
+            const dictionary& solverDict
+        );
+
+
+    // Destructor
+
+        virtual ~TPCG()
+        {}
+
+
+    // Member Functions
+
+        //- Solve the matrix with this solver
+        virtual typename LduMatrix<Type, DType, LUType>::solverPerformance solve
+        (
+            Field<Type>& psi
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "TPCG.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.C

@@ -0,0 +1,153 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "SmoothSolver.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+Foam::SmoothSolver<Type, DType, LUType>::SmoothSolver
+(
+    const word& fieldName,
+    const LduMatrix<Type, DType, LUType>& matrix,
+    const dictionary& solverDict
+)
+:
+    LduMatrix<Type, DType, LUType>::solver
+    (
+        fieldName,
+        matrix,
+        solverDict
+    ),
+    nSweeps_(1)
+{
+    readControls();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type, class DType, class LUType>
+void Foam::SmoothSolver<Type, DType, LUType>::readControls()
+{
+    LduMatrix<Type, DType, LUType>::solver::readControls();
+    readControl(this->controlDict_, nSweeps_, "nSweeps");
+}
+
+
+template<class Type, class DType, class LUType>
+typename Foam::LduMatrix<Type, DType, LUType>::solverPerformance
+Foam::SmoothSolver<Type, DType, LUType>::solve(Field<Type>& psi) const
+{
+    // --- Setup class containing solver performance data
+    typename LduMatrix<Type, DType, LUType>::solverPerformance solverPerf
+    (
+        typeName,
+        this->fieldName_
+    );
+
+    // If the nSweeps_ is negative do a fixed number of sweeps
+    if (nSweeps_ < 0)
+    {
+        autoPtr<typename LduMatrix<Type, DType, LUType>::smoother>
+        smootherPtr = LduMatrix<Type, DType, LUType>::smoother::New
+        (
+            this->fieldName_,
+            this->matrix_,
+            this->controlDict_
+        );
+
+        smootherPtr->smooth(psi, -nSweeps_);
+
+        solverPerf.nIterations() -= nSweeps_;
+    }
+    else
+    {
+        Type normFactor = pTraits<Type>::zero;
+
+        {
+            Field<Type> Apsi(psi.size());
+            Field<Type> temp(psi.size());
+
+            // Calculate A.psi
+            this->matrix_.Amul(Apsi, psi);
+
+            // Calculate normalisation factor
+            normFactor = this->normFactor(psi, Apsi, temp);
+
+            // Calculate residual magnitude
+            solverPerf.initialResidual() = cmptDivide
+            (
+                gSumCmptMag(this->matrix_.source() - Apsi),
+                normFactor
+            );
+            solverPerf.finalResidual() = solverPerf.initialResidual();
+        }
+
+        if (LduMatrix<Type, DType, LUType>::debug >= 2)
+        {
+            Info<< "   Normalisation factor = " << normFactor << endl;
+        }
+
+
+        // Check convergence, solve if not converged
+        if (!solverPerf.converged(this->tolerance_, this->relTol_))
+        {
+            autoPtr<typename LduMatrix<Type, DType, LUType>::smoother>
+            smootherPtr = LduMatrix<Type, DType, LUType>::smoother::New
+            (
+                this->fieldName_,
+                this->matrix_,
+                this->controlDict_
+            );
+
+            // Smoothing loop
+            do
+            {
+                smootherPtr->smooth
+                (
+                    psi,
+                    nSweeps_
+                );
+
+                // Calculate the residual to check convergence
+                solverPerf.finalResidual() = cmptDivide
+                (
+                    gSumCmptMag(this->matrix_.residual(psi)),
+                    normFactor
+                );
+            } while
+            (
+                (solverPerf.nIterations() += nSweeps_) < this->maxIter_
+             && !(solverPerf.converged(this->tolerance_, this->relTol_))
+            );
+        }
+    }
+
+    return solverPerf;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.H

@@ -0,0 +1,110 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::SmoothSolver
+
+Description
+    Iterative solver for symmetric and assymetric matrices which uses a
+    run-time selected smoother e.g. GaussSeidel to converge the solution to
+    the required tolerance.  To improve efficiency, the residual is evaluated
+    after every nSweeps smoothing iterations.
+
+SourceFiles
+    SmoothSolver.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef SmoothSolver_H
+#define SmoothSolver_H
+
+#include "lduMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class SmoothSolver Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type, class DType, class LUType>
+class SmoothSolver
+:
+    public LduMatrix<Type, DType, LUType>::solver
+{
+
+protected:
+
+    // Protected data
+
+        //- Number of sweeps before the evaluation of residual
+        label nSweeps_;
+
+        //- Read the control parameters from the controlDict_
+        virtual void readControls();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("SmoothSolver");
+
+
+    // Constructors
+
+        //- Construct from matrix components and solver data dictionary
+        SmoothSolver
+        (
+            const word& fieldName,
+            const LduMatrix<Type, DType, LUType>& matrix,
+            const dictionary& solverDict
+        );
+
+
+    // Member Functions
+
+        //- Solve the matrix with this solver
+        virtual typename LduMatrix<Type, DType, LUType>::solverPerformance solve
+        (
+            Field<Type>& psi
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "SmoothSolver.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/matrices/LduMatrix/Solvers/lduSolvers.C

@@ -0,0 +1,57 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "TPCG.H"
+#include "TPBiCG.H"
+#include "SmoothSolver.H"
+#include "fieldTypes.H"
+
+#define makeLduSolvers(Type, DType, LUType)                                   \
+                                                                              \
+    makeLduSolver(DiagonalSolver, Type, DType, LUType);                       \
+    makeLduSymSolver(DiagonalSolver, Type, DType, LUType);                    \
+    makeLduAsymSolver(DiagonalSolver, Type, DType, LUType);                   \
+                                                                              \
+    makeLduSolver(TPCG, Type, DType, LUType);                                 \
+    makeLduSymSolver(TPCG, Type, DType, LUType);                              \
+                                                                              \
+    makeLduSolver(TPBiCG, Type, DType, LUType);                               \
+    makeLduAsymSolver(TPBiCG, Type, DType, LUType);                           \
+                                                                              \
+    makeLduSolver(SmoothSolver, Type, DType, LUType);                         \
+    makeLduSymSolver(SmoothSolver, Type, DType, LUType);                      \
+    makeLduAsymSolver(SmoothSolver, Type, DType, LUType);
+
+namespace Foam
+{
+    makeLduSolvers(scalar, scalar, scalar);
+    makeLduSolvers(vector, scalar, scalar);
+    makeLduSolvers(sphericalTensor, scalar, scalar);
+    makeLduSolvers(symmTensor, scalar, scalar);
+    makeLduSolvers(tensor, scalar, scalar);
+};
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterfaceFields/lduInterfaceField/lduInterfaceField.H

@@ -115,9 +115,8 @@ public:
             //- Initialise neighbour matrix update
             virtual void initInterfaceMatrixUpdate
             (
-                const scalarField&,
                 scalarField&,
-                const lduMatrix&,
+                const scalarField&,
                 const scalarField&,
                 const direction,
                 const Pstream::commsTypes commsType
@@ -145,9 +144,8 @@ public:
             //- Update result field based on interface functionality
             virtual void updateInterfaceMatrix
             (
-                const scalarField&,
                 scalarField&,
-                const lduMatrix&,
+                const scalarField&,
                 const scalarField&,
                 const direction,
                 const Pstream::commsTypes commsType

src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixUpdateMatrixInterfaces.C

@@ -48,9 +48,8 @@ void Foam::lduMatrix::initMatrixInterfaces
             {
                 interfaces[interfaceI].initInterfaceMatrixUpdate
                 (
-                    psiif,
                     result,
-                    *this,
+                    psiif,
                     coupleCoeffs[interfaceI],
                     cmpt,
                     Pstream::defaultCommsType
@@ -75,9 +74,8 @@ void Foam::lduMatrix::initMatrixInterfaces
             {
                 interfaces[interfaceI].initInterfaceMatrixUpdate
                 (
-                    psiif,
                     result,
-                    *this,
+                    psiif,
                     coupleCoeffs[interfaceI],
                     cmpt,
                     Pstream::blocking
@@ -112,9 +110,8 @@ void Foam::lduMatrix::updateMatrixInterfaces
             {
                 interfaces[interfaceI].updateInterfaceMatrix
                 (
-                    psiif,
                     result,
-                    *this,
+                    psiif,
                     coupleCoeffs[interfaceI],
                     cmpt,
                     Pstream::defaultCommsType
@@ -141,9 +138,8 @@ void Foam::lduMatrix::updateMatrixInterfaces
                         {
                             interfaces[interfaceI].updateInterfaceMatrix
                             (
-                                psiif,
                                 result,
-                                *this,
+                                psiif,
                                 coupleCoeffs[interfaceI],
                                 cmpt,
                                 Pstream::defaultCommsType
@@ -193,9 +189,8 @@ void Foam::lduMatrix::updateMatrixInterfaces
             {
                 interfaces[interfaceI].updateInterfaceMatrix
                 (
-                    psiif,
                     result,
-                    *this,
+                    psiif,
                     coupleCoeffs[interfaceI],
                     cmpt,
                     Pstream::defaultCommsType
@@ -218,9 +213,8 @@ void Foam::lduMatrix::updateMatrixInterfaces
                 {
                     interfaces[interfaceI].initInterfaceMatrixUpdate
                     (
-                        psiif,
                         result,
-                        *this,
+                        psiif,
                         coupleCoeffs[interfaceI],
                         cmpt,
                         Pstream::scheduled
@@ -230,9 +224,8 @@ void Foam::lduMatrix::updateMatrixInterfaces
                 {
                     interfaces[interfaceI].updateInterfaceMatrix
                     (
-                        psiif,
                         result,
-                        *this,
+                        psiif,
                         coupleCoeffs[interfaceI],
                         cmpt,
                         Pstream::scheduled
@@ -254,9 +247,8 @@ void Foam::lduMatrix::updateMatrixInterfaces
             {
                 interfaces[interfaceI].updateInterfaceMatrix
                 (
-                    psiif,
                     result,
-                    *this,
+                    psiif,
                     coupleCoeffs[interfaceI],
                     cmpt,
                     Pstream::blocking

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/cyclicGAMGInterfaceField/cyclicGAMGInterfaceField.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -72,9 +72,8 @@ Foam::cyclicGAMGInterfaceField::~cyclicGAMGInterfaceField()
 
 void Foam::cyclicGAMGInterfaceField::updateInterfaceMatrix
 (
-    const scalarField& psiInternal,
     scalarField& result,
-    const lduMatrix&,
+    const scalarField& psiInternal,
     const scalarField& coeffs,
     const direction cmpt,
     const Pstream::commsTypes

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/cyclicGAMGInterfaceField/cyclicGAMGInterfaceField.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -110,9 +110,8 @@ public:
             //- Update result field based on interface functionality
             virtual void updateInterfaceMatrix
             (
-                const scalarField& psiInternal,
                 scalarField& result,
-                const lduMatrix&,
+                const scalarField& psiInternal,
                 const scalarField& coeffs,
                 const direction cmpt,
                 const Pstream::commsTypes commsType

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/processorCyclicGAMGInterfaceField/processorCyclicGAMGInterfaceField.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -63,9 +63,8 @@ Foam::processorCyclicGAMGInterfaceField::~processorCyclicGAMGInterfaceField()
 
 //void Foam::processorCyclicGAMGInterfaceField::initInterfaceMatrixUpdate
 //(
-//    const scalarField& psiInternal,
 //    scalarField&,
-//    const lduMatrix&,
+//    const scalarField&,
 //    const scalarField&,
 //    const direction,
 //    const Pstream::commsTypes commsType
@@ -81,9 +80,8 @@ Foam::processorCyclicGAMGInterfaceField::~processorCyclicGAMGInterfaceField()
 //
 //void Foam::processorCyclicGAMGInterfaceField::updateInterfaceMatrix
 //(
-//    const scalarField&,
 //    scalarField& result,
-//    const lduMatrix&,
+//    const scalarField&,
 //    const scalarField& coeffs,
 //    const direction cmpt,
 //    const Pstream::commsTypes commsType

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/processorCyclicGAMGInterfaceField/processorCyclicGAMGInterfaceField.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -111,9 +111,8 @@ public:
 //            //- Initialise neighbour matrix update
 //            virtual void initInterfaceMatrixUpdate
 //            (
-//                const scalarField& psiInternal,
 //                scalarField& result,
-//                const lduMatrix& m,
+//                const scalarField& psiInternal,
 //                const scalarField& coeffs,
 //                const direction cmpt,
 //                const Pstream::commsTypes commsType
@@ -122,9 +121,8 @@ public:
 //            //- Update result field based on interface functionality
 //            virtual void updateInterfaceMatrix
 //            (
-//                const scalarField& psiInternal,
 //                scalarField& result,
-//                const lduMatrix&,
+//                const scalarField& psiInternal,
 //                const scalarField& coeffs,
 //                const direction cmpt,
 //                const Pstream::commsTypes commsType

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/processorGAMGInterfaceField/processorGAMGInterfaceField.C

@@ -72,9 +72,8 @@ Foam::processorGAMGInterfaceField::~processorGAMGInterfaceField()
 
 void Foam::processorGAMGInterfaceField::initInterfaceMatrixUpdate
 (
-    const scalarField& psiInternal,
     scalarField&,
-    const lduMatrix&,
+    const scalarField& psiInternal,
     const scalarField&,
     const direction,
     const Pstream::commsTypes commsType
@@ -117,9 +116,8 @@ void Foam::processorGAMGInterfaceField::initInterfaceMatrixUpdate
 
 void Foam::processorGAMGInterfaceField::updateInterfaceMatrix
 (
-    const scalarField&,
     scalarField& result,
-    const lduMatrix&,
+    const scalarField&,
     const scalarField& coeffs,
     const direction cmpt,
     const Pstream::commsTypes commsType

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/processorGAMGInterfaceField/processorGAMGInterfaceField.H

@@ -126,9 +126,8 @@ public:
             //- Initialise neighbour matrix update
             virtual void initInterfaceMatrixUpdate
             (
-                const scalarField& psiInternal,
                 scalarField& result,
-                const lduMatrix& m,
+                const scalarField& psiInternal,
                 const scalarField& coeffs,
                 const direction cmpt,
                 const Pstream::commsTypes commsType
@@ -137,9 +136,8 @@ public:
             //- Update result field based on interface functionality
             virtual void updateInterfaceMatrix
             (
-                const scalarField& psiInternal,
                 scalarField& result,
-                const lduMatrix&,
+                const scalarField& psiInternal,
                 const scalarField& coeffs,
                 const direction cmpt,
                 const Pstream::commsTypes commsType

src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPoints.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2012 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/dynamicMesh/extrudePatchMesh/extrudePatchMesh.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2012 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/dynamicMesh/extrudePatchMesh/extrudePatchMesh.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2012 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/finiteVolume/fields/fvPatchFields/basic/coupled/coupledFvPatchField.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -178,9 +178,8 @@ public:
             //- Update result field based on interface functionality
             virtual void updateInterfaceMatrix
             (
-                const scalarField& psiInternal,
                 scalarField& result,
-                const lduMatrix&,
+                const scalarField& psiInternal,
                 const scalarField& coeffs,
                 const direction,
                 const Pstream::commsTypes commsType

src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -191,9 +191,8 @@ const
 template<class Type>
 void cyclicFvPatchField<Type>::updateInterfaceMatrix
 (
-    const scalarField& psiInternal,
     scalarField& result,
-    const lduMatrix&,
+    const scalarField& psiInternal,
     const scalarField& coeffs,
     const direction cmpt,
     const Pstream::commsTypes

src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -159,9 +159,8 @@ public:
             //- Update result field based on interface functionality
             virtual void updateInterfaceMatrix
             (
-                const scalarField& psiInternal,
                 scalarField& result,
-                const lduMatrix&,
+                const scalarField& psiInternal,
                 const scalarField& coeffs,
                 const direction cmpt,
                 const Pstream::commsTypes commsType

src/finiteVolume/fields/fvPatchFields/constraint/cyclicAMI/cyclicAMIFvPatchField.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -202,9 +202,8 @@ Foam::cyclicAMIFvPatchField<Type>::neighbourPatchField() const
 template<class Type>
 void Foam::cyclicAMIFvPatchField<Type>::updateInterfaceMatrix
 (
-    const scalarField& psiInternal,
     scalarField& result,
-    const lduMatrix&,
+    const scalarField& psiInternal,
     const scalarField& coeffs,
     const direction cmpt,
     const Pstream::commsTypes

src/finiteVolume/fields/fvPatchFields/constraint/cyclicAMI/cyclicAMIFvPatchField.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -160,9 +160,8 @@ public:
             //- Update result field based on interface functionality
             virtual void updateInterfaceMatrix
             (
-                const scalarField& psiInternal,
                 scalarField& result,
-                const lduMatrix&,
+                const scalarField& psiInternal,
                 const scalarField& coeffs,
                 const direction cmpt,
                 const Pstream::commsTypes commsType

src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -137,9 +137,8 @@ tmp<Field<Type> > jumpCyclicFvPatchField<Type>::patchNeighbourField() const
 template<class Type>
 void jumpCyclicFvPatchField<Type>::updateInterfaceMatrix
 (
-    const scalarField& psiInternal,
     scalarField& result,
-    const lduMatrix&,
+    const scalarField& psiInternal,
     const scalarField& coeffs,
     const direction cmpt,
     const Pstream::commsTypes

src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -120,9 +120,8 @@ public:
             //- Update result field based on interface functionality
             virtual void updateInterfaceMatrix
             (
-                const scalarField& psiInternal,
                 scalarField& result,
-                const lduMatrix&,
+                const scalarField& psiInternal,
                 const scalarField& coeffs,
                 const direction cmpt,
                 const Pstream::commsTypes commsType

src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchField.C

@@ -311,9 +311,8 @@ tmp<Field<Type> > processorFvPatchField<Type>::snGrad() const
 template<class Type>
 void processorFvPatchField<Type>::initInterfaceMatrixUpdate
 (
-    const scalarField& psiInternal,
     scalarField&,
-    const lduMatrix&,
+    const scalarField& psiInternal,
     const scalarField&,
     const direction,
     const Pstream::commsTypes commsType
@@ -368,9 +367,8 @@ void processorFvPatchField<Type>::initInterfaceMatrixUpdate
 template<class Type>
 void processorFvPatchField<Type>::updateInterfaceMatrix
 (
-    const scalarField&,
     scalarField& result,
-    const lduMatrix&,
+    const scalarField&,
     const scalarField& coeffs,
     const direction cmpt,
     const Pstream::commsTypes commsType

src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchField.H

@@ -187,9 +187,8 @@ public:
             //- Initialise neighbour matrix update
             virtual void initInterfaceMatrixUpdate
             (
-                const scalarField& psiInternal,
                 scalarField& result,
-                const lduMatrix& m,
+                const scalarField& psiInternal,
                 const scalarField& coeffs,
                 const direction cmpt,
                 const Pstream::commsTypes commsType
@@ -201,9 +200,8 @@ public:
             //- Update result field based on interface functionality
             virtual void updateInterfaceMatrix
             (
-                const scalarField& psiInternal,
                 scalarField& result,
-                const lduMatrix& m,
+                const scalarField& psiInternal,
                 const scalarField& coeffs,
                 const direction cmpt,
                 const Pstream::commsTypes commsType

src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchScalarField.C

@@ -35,9 +35,8 @@ namespace Foam
 template<>
 void processorFvPatchField<scalar>::initInterfaceMatrixUpdate
 (
-    const scalarField& psiInternal,
     scalarField&,
-    const lduMatrix&,
+    const scalarField& psiInternal,
     const scalarField&,
     const direction,
     const Pstream::commsTypes commsType
@@ -92,9 +91,8 @@ void processorFvPatchField<scalar>::initInterfaceMatrixUpdate
 template<>
 void processorFvPatchField<scalar>::updateInterfaceMatrix
 (
-    const scalarField&,
     scalarField& result,
-    const lduMatrix&,
+    const scalarField&,
     const scalarField& coeffs,
     const direction,
     const Pstream::commsTypes commsType

src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchScalarField.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,9 +38,8 @@ namespace Foam
 template<>
 void processorFvPatchField<scalar>::initInterfaceMatrixUpdate
 (
-    const scalarField&,
     scalarField&,
-    const lduMatrix&,
+    const scalarField&,
     const scalarField&,
     const direction,
     const Pstream::commsTypes commsType
@@ -50,9 +49,8 @@ void processorFvPatchField<scalar>::initInterfaceMatrixUpdate
 template<>
 void processorFvPatchField<scalar>::updateInterfaceMatrix
 (
-    const scalarField&,
     scalarField& result,
-    const lduMatrix&,
+    const scalarField&,
     const scalarField& coeffs,
     const direction,
     const Pstream::commsTypes commsType

src/finiteVolume/fields/fvPatchFields/constraint/processorCyclic/processorCyclicFvPatchField.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -227,9 +227,8 @@ processorCyclicFvPatchField<Type>::~processorCyclicFvPatchField()
 //template<class Type>
 //void processorCyclicFvPatchField<Type>::initInterfaceMatrixUpdate
 //(
-//    const scalarField& psiInternal,
 //    scalarField&,
-//    const lduMatrix&,
+//    const scalarField& psiInternal,
 //    const scalarField&,
 //    const direction,
 //    const Pstream::commsTypes commsType
@@ -246,9 +245,8 @@ processorCyclicFvPatchField<Type>::~processorCyclicFvPatchField()
 //template<class Type>
 //void processorCyclicFvPatchField<Type>::updateInterfaceMatrix
 //(
-//    const scalarField&,
 //    scalarField& result,
-//    const lduMatrix&,
+//    const scalarField&,
 //    const scalarField& coeffs,
 //    const direction cmpt,
 //    const Pstream::commsTypes commsType

src/finiteVolume/fields/fvPatchFields/constraint/processorCyclic/processorCyclicFvPatchScalarField.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,9 +38,8 @@ namespace Foam
 //template<>
 //void processorCyclicFvPatchField<scalar>::initInterfaceMatrixUpdate
 //(
-//    const scalarField&,
 //    scalarField&,
-//    const lduMatrix&,
+//    const scalarField&,
 //    const scalarField&,
 //    const direction,
 //    const Pstream::commsTypes commsType
@@ -50,9 +49,8 @@ namespace Foam
 //template<>
 //void processorCyclicFvPatchField<scalar>::updateInterfaceMatrix
 //(
-//    const scalarField&,
 //    scalarField& result,
-//    const lduMatrix&,
+//    const scalarField&,
 //    const scalarField& coeffs,
 //    const direction,
 //    const Pstream::commsTypes commsType

src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -142,7 +142,8 @@ class fvMatrix
         mutable GeometricField<Type, fvsPatchField, surfaceMesh>
             *faceFluxCorrectionPtr_;
 
-protected:
+// ***HGW for testing LduMatrix protected:
+public:
 
     //- Declare friendship with the fvSolver class
     friend class fvSolver;

src/meshTools/AMIInterpolation/GAMG/interfaceFields/cyclicAMIGAMGInterfaceField/cyclicAMIGAMGInterfaceField.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -72,9 +72,8 @@ Foam::cyclicAMIGAMGInterfaceField::~cyclicAMIGAMGInterfaceField()
 
 void Foam::cyclicAMIGAMGInterfaceField::updateInterfaceMatrix
 (
-    const scalarField& psiInternal,
     scalarField& result,
-    const lduMatrix&,
+    const scalarField& psiInternal,
     const scalarField& coeffs,
     const direction cmpt,
     const Pstream::commsTypes

src/meshTools/AMIInterpolation/GAMG/interfaceFields/cyclicAMIGAMGInterfaceField/cyclicAMIGAMGInterfaceField.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-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -110,9 +110,8 @@ public:
             //- Update result field based on interface functionality
             virtual void updateInterfaceMatrix
             (
-                const scalarField& psiInternal,
                 scalarField& result,
-                const lduMatrix&,
+                const scalarField& psiInternal,
                 const scalarField& coeffs,
                 const direction cmpt,
                 const Pstream::commsTypes commsType

src/triSurface/triSurface/triSurface.C

@@ -386,7 +386,7 @@ bool Foam::triSurface::read
     }
     else if (ext == "stlb")
     {
-        return readSTL(name);
+        return readSTLBINARY(name);
     }
     else if (ext == "gts")
     {