Feature: Comulative release fixes. Author: Hrvoje Jasak. Merge: Hrvoje Jasak.

applications/utilities/surface/surfaceAutoPatch/Make/options

@@ -1,5 +1,5 @@
 EXE_INC = \
-    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
 
 EXE_LIBS = \
-    -lmeshTools \
+    -lmeshTools

applications/utilities/surface/surfaceClean/Make/options

@@ -2,4 +2,4 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude
 
 EXE_LIBS = \
-    -lmeshTools \
+    -lmeshTools

applications/utilities/surface/surfaceFeatureConvert/Make/options

@@ -4,4 +4,4 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lmeshTools \
-    -ledgeMesh \
+    -ledgeMesh

applications/utilities/surface/surfaceFeatureExtract/Make/options

@@ -1,6 +1,6 @@
 EXE_INC = \
     -I$(LIB_SRC)/cfdTools/general/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
 
 EXE_LIBS = \
-    -lmeshTools \
+    -lmeshTools

applications/utilities/surface/surfaceMeshConvert/Make/options

@@ -2,4 +2,6 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude
 
-EXE_LIBS = -lmeshTools -lsurfMesh
+EXE_LIBS = \
+    -lmeshTools \
+    -lsurfMesh

applications/utilities/surface/surfaceMeshExport/Make/options

@@ -2,4 +2,6 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude
 
-EXE_LIBS = -lmeshTools -lsurfMesh
+EXE_LIBS = \
+    -lmeshTools \
+    -lsurfMesh

applications/utilities/surface/surfaceMeshImport/Make/options

@@ -2,4 +2,6 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude
 
-EXE_LIBS = -lmeshTools -lsurfMesh
+EXE_LIBS = \
+    -lmeshTools \
+    -lsurfMesh

applications/utilities/surface/surfaceMeshTriangulate/Make/options

@@ -1,6 +1,5 @@
 EXE_INC = \
-    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
 
 EXE_LIBS = \
-    -lmeshTools \
-
+    -lmeshTools

applications/utilities/surface/surfaceOrient/Make/options

@@ -2,4 +2,4 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude
 
 EXE_LIBS = \
-    -lmeshTools \
+    -lmeshTools

applications/utilities/surface/surfacePointMerge/Make/options

@@ -2,4 +2,4 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
 
 EXE_LIBS = \
-    -lmeshTools \
+    -lmeshTools

applications/utilities/surface/surfaceRedistributePar/Make/options

@@ -2,4 +2,4 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude
 
 EXE_LIBS = \
-    -lmeshTools \
+    -lmeshTools

applications/utilities/surface/surfaceRefineRedGreen/Make/options

@@ -1,5 +1,5 @@
 EXE_INC = \
-    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
 
 EXE_LIBS = \
-    -lmeshTools \
+    -lmeshTools

applications/utilities/surface/surfaceSplitByPatch/Make/options

@@ -3,4 +3,3 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lmeshTools
-

applications/utilities/surface/surfaceSplitNonManifolds/Make/options

@@ -1,5 +1,5 @@
 EXE_INC = \
-    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
 
 EXE_LIBS = \
-    -lmeshTools \
+    -lmeshTools

src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C

@@ -68,8 +68,8 @@ reconstruct
                 IOobject::NO_READ,
                 IOobject::NO_WRITE
             ),
-            inv(surfaceSum(sqr(mesh.Sf())/mesh.magSf()))
-          & surfaceSum((mesh.Sf()/mesh.magSf())*ssf),
+            mesh,
+            ssf.dimensions()/dimArea,
             zeroGradientFvPatchField<GradType>::typeName
         )
     );
@@ -87,9 +87,14 @@ reconstruct
     GeometricField<GradType, fvPatchField, volMesh> fluxTimesNormal =
         surfaceSum((mesh.Sf()/mesh.magSf())*ssf);
 
+    // Note: hinv inverse must be used to stabilise the inverse on bad meshes
+    // HJ, 19/Aug/2015
     reconField.internalField() =
     (
-        inv(surfaceSum(sqr(mesh.Sf())/mesh.magSf())().internalField())
+        hinv
+        (
+            surfaceSum(sqr(mesh.Sf())/mesh.magSf())().internalField()
+        )
       & fluxTimesNormal.internalField()
     );
 

src/foam/Make/files

@@ -689,6 +689,11 @@ matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/blockGaussSeidelP
 matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/scalarBlockCholeskyPrecon.C
 matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/tensorBlockCholeskyPrecon.C
 matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/blockCholeskyPrecons.C
+
+matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/scalarBlockDiagCholeskyPrecon.C
+matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/tensorBlockDiagCholeskyPrecon.C
+matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/blockDiagCholeskyPrecons.C
+
 matrices/blockLduMatrix/BlockLduPrecons/BlockAmgPrecon/blockAmgPrecons.C
 
 matrices/blockLduMatrix/BlockLduSmoothers/BlockLduSmoother/blockLduSmoothers.C

src/foam/fields/CoeffField/DecoupledCoeffFieldFunctions.H

@@ -37,7 +37,10 @@ namespace Foam
 /* * * * * * * * * * * * * * * * Global functions  * * * * * * * * * * * * * */
 
 template<class Type>
-tmp<DecoupledCoeffField<Type> > inv(const DecoupledCoeffField<Type>& f);
+tmp<DecoupledCoeffField<Type> > inv
+(
+    const DecoupledCoeffField<Type>& f
+);
 
 
 template<class Type>

src/foam/interpolations/GGIInterpolation/GGIInterpolationQuickRejectTests.C

@@ -181,7 +181,7 @@ void GGIInterpolation<MasterPatch, SlavePatch>::findNeighbours3D
 
         scalar tmpValue = Foam::magSqr(bbSlave.max() - bbSlave.min())/4.0;
 
-        // We will compare squared distances, so save the sqrt() if value > 1.0.
+        // We compare squared distances, so save the sqrt() if value > 1.0.
         if (tmpValue < 1.0)
         {
             // Take the sqrt, otherwise, we underestimate the radius
@@ -566,7 +566,6 @@ void GGIInterpolation<MasterPatch, SlavePatch>::findNeighboursBBOctree
 }
 
 
-
 // Projects a list of points onto a plane located at planeOrig,
 // oriented along planeNormal.  Return the projected points in a
 // pointField, and the normal distance of each points from the

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/BlockCholeskyPrecon.C

@@ -220,6 +220,9 @@ void Foam::BlockCholeskyPrecon<Type>::calcPreconDiag()
             }
         }
     }
+
+    // Invert the diagonal
+    preconDiag_ = inv(preconDiag_);
 }
 
 
@@ -250,12 +253,6 @@ void Foam::BlockCholeskyPrecon<Type>::diagMultiply
                 upper[coeffI]
             );
     }
-
-    // Invert the diagonal for future use
-    forAll (dDiag, i)
-    {
-        dDiag[i] = mult.inverse(dDiag[i]);
-    }
 }
 
 
@@ -286,12 +283,6 @@ void Foam::BlockCholeskyPrecon<Type>::diagMultiplyCoeffT
                 upper[coeffI]
             );
     }
-
-    // Invert the diagonal for future use
-    forAll (dDiag, i)
-    {
-        dDiag[i] = mult.inverse(dDiag[i]);
-    }
 }
 
 
@@ -323,12 +314,6 @@ void Foam::BlockCholeskyPrecon<Type>::diagMultiply
                 upper[coeffI]
             );
     }
-
-    // Invert the diagonal for future use
-    forAll (dDiag, i)
-    {
-        dDiag[i] = mult.inverse(dDiag[i]);
-    }
 }
 
 

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/BlockCholeskyPreconDecoupled.C

@@ -126,6 +126,11 @@ void Foam::BlockCholeskyPrecon<Type>::calcDecoupledPreconDiag()
             }
         }
     }
+
+    // Invert the diagonal
+    // Note: since square diag type does not exist, simple inversion
+    // covers all cases
+    preconDiag_ = inv(preconDiag_);
 }
 
 

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/BlockDiagCholeskyPrecon.H

@@ -0,0 +1,235 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     3.2
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    BlockDiagCholeskyPrecon
+
+Description
+    Incomplete Cholesky preconditioning with no fill-in,
+    using the diagonal-of-diagonal for ILU decomposition.
+    Currently broken
+
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved.
+
+SourceFiles
+    BlockDiagCholeskyPrecon.C
+    BlockDiagCholeskyPreconDecoupled.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef BlockDiagCholeskyPrecon_H
+#define BlockDiagCholeskyPrecon_H
+
+#include "BlockLduPrecon.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class BlockDiagCholeskyPrecon Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class BlockDiagCholeskyPrecon
+:
+    public BlockLduPrecon<Type>
+{
+    // Private Data
+
+        //- Preconditioned diagonal
+        mutable CoeffField<Type> preconDiag_;
+
+        //- Off-diag part of diagonal
+        CoeffField<Type> LUDiag_;
+
+        //- Temporary space for updated decoupled source
+        //  Initialised with zero size and resized on first use
+        mutable Field<Type> bPlusLU_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        BlockDiagCholeskyPrecon(const BlockDiagCholeskyPrecon&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const BlockDiagCholeskyPrecon&);
+
+        //- Precondition the diagonal
+        void calcPreconDiag();
+
+        // Diagonal multiplication, symmetric matrix
+        template<class DiagType, class ULType>
+        void diagMultiply
+        (
+            Field<DiagType>& dDiag,
+            const Field<ULType>& upper
+        );
+
+        //- Diagonal multiplication with transpose upper square coeff
+        //  for the symmetric matrix
+        template<class DiagType, class ULType>
+        void diagMultiplyCoeffT
+        (
+            Field<DiagType>& dDiag,
+            const Field<ULType>& upper
+        );
+
+        //- Diagonal multiplication, asymmetric matrix
+        template<class DiagType, class ULType>
+        void diagMultiply
+        (
+            Field<DiagType>& dDiag,
+            const Field<ULType>& lower,
+            const Field<ULType>& upper
+        );
+
+        //- ILU multiplication, symmetric matrix
+        template<class DiagType, class ULType>
+        void ILUmultiply
+        (
+            Field<Type>& x,
+            const Field<DiagType>& dDiag,
+            const Field<ULType>& upper,
+            const Field<Type>& b
+        ) const;
+
+        //- ILU multiplication, with transpose upper square coeff
+        //  for a symmetric matrix
+        template<class DiagType, class ULType>
+        void ILUmultiplyCoeffT
+        (
+            Field<Type>& x,
+            const Field<DiagType>& dDiag,
+            const Field<ULType>& upper,
+            const Field<Type>& b
+        ) const;
+
+        //- ILU multiplication, asymmetric matrix
+        template<class DiagType, class ULType>
+        void ILUmultiply
+        (
+            Field<Type>& x,
+            const Field<DiagType>& dDiag,
+            const Field<ULType>& lower,
+            const Field<ULType>& upper,
+            const Field<Type>& b
+        ) const;
+
+        //- ILU multiplication transposed asymmetric matrix
+        template<class DiagType, class ULType>
+        void ILUmultiplyTranspose
+        (
+            Field<Type>& x,
+            const Field<DiagType>& dDiag,
+            const Field<ULType>& lower,
+            const Field<ULType>& upper,
+            const Field<Type>& b
+        ) const;
+
+
+        // Decoupled operations, used in template specialisation
+
+            //- Precondition the diagonal, decoupled version
+            void calcDecoupledPreconDiag();
+
+            //- Execute preconditioning, decoupled version
+            void decoupledPrecondition
+            (
+                Field<Type>& x,
+                const Field<Type>& b
+            ) const;
+
+            //- Execute preconditioning with matrix transpose,
+            //  decoupled version
+            void decoupledPreconditionT
+            (
+                Field<Type>& xT,
+                const Field<Type>& bT
+            ) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("Cholesky");
+
+
+    // Constructors
+
+        //- Construct from matrix for smoother use
+        BlockDiagCholeskyPrecon
+        (
+            const BlockLduMatrix<Type>& matrix
+        );
+
+        //- Construct from components
+        BlockDiagCholeskyPrecon
+        (
+            const BlockLduMatrix<Type>& matrix,
+            const dictionary& dict
+        );
+
+
+    // Destructor
+
+        virtual ~BlockDiagCholeskyPrecon();
+
+
+    // Member Functions
+
+        //- Execute preconditioning
+        virtual void precondition
+        (
+            Field<Type>& x,
+            const Field<Type>& b
+        ) const;
+
+        //- Execute preconditioning with matrix transpose
+        virtual void preconditionT
+        (
+            Field<Type>& xT,
+            const Field<Type>& bT
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "BlockDiagCholeskyPrecon.C"
+#   include "BlockDiagCholeskyPreconDecoupled.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/BlockDiagCholeskyPreconDecoupled.C

@@ -0,0 +1,332 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     3.2
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+#include "BlockDiagCholeskyPrecon.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+void Foam::BlockDiagCholeskyPrecon<Type>::calcDecoupledPreconDiag()
+{
+    typedef CoeffField<Type> TypeCoeffField;
+
+    // Note: Assuming lower and upper triangle have the same active type
+
+    if (this->matrix_.symmetric())
+    {
+        const TypeCoeffField& UpperCoeff = this->matrix_.upper();
+
+        if (preconDiag_.activeType() == blockCoeffBase::SCALAR)
+        {
+            if (UpperCoeff.activeType() == blockCoeffBase::SCALAR)
+            {
+                diagMultiply
+                (
+                    preconDiag_.asScalar(),
+                    UpperCoeff.asScalar()
+                );
+            }
+            else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR)
+            {
+                diagMultiply
+                (
+                    preconDiag_.asLinear(),
+                    UpperCoeff.asLinear()
+                );
+            }
+        }
+        else if (preconDiag_.activeType() == blockCoeffBase::LINEAR)
+        {
+            if (UpperCoeff.activeType() == blockCoeffBase::SCALAR)
+            {
+                diagMultiply
+                (
+                    preconDiag_.asLinear(),
+                    UpperCoeff.asScalar()
+                );
+            }
+            else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR)
+            {
+                diagMultiply
+                (
+                    preconDiag_.asLinear(),
+                    UpperCoeff.asLinear()
+                );
+            }
+        }
+    }
+    else // Asymmetric matrix
+    {
+        const TypeCoeffField& LowerCoeff = this->matrix_.lower();
+        const TypeCoeffField& UpperCoeff = this->matrix_.upper();
+
+        if (preconDiag_.activeType() == blockCoeffBase::SCALAR)
+        {
+            if (UpperCoeff.activeType() == blockCoeffBase::SCALAR)
+            {
+                diagMultiply
+                (
+                    preconDiag_.asScalar(),
+                    LowerCoeff.asScalar(),
+                    UpperCoeff.asScalar()
+                );
+            }
+            else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR)
+            {
+                diagMultiply
+                (
+                    preconDiag_.asLinear(),
+                    LowerCoeff.asLinear(),
+                    UpperCoeff.asLinear()
+                );
+            }
+        }
+        else if (preconDiag_.activeType() == blockCoeffBase::LINEAR)
+        {
+            if (UpperCoeff.activeType() == blockCoeffBase::SCALAR)
+            {
+                diagMultiply
+                (
+                    preconDiag_.asLinear(),
+                    LowerCoeff.asScalar(),
+                    UpperCoeff.asScalar()
+                );
+            }
+            else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR)
+            {
+                diagMultiply
+                (
+                    preconDiag_.asLinear(),
+                    LowerCoeff.asLinear(),
+                    UpperCoeff.asLinear()
+                );
+            }
+        }
+    }
+
+    // Invert the diagonal
+    // Note: since square diag type does not exist, simple inversion
+    // covers all cases
+    preconDiag_ = inv(preconDiag_);
+}
+
+
+template<class Type>
+void Foam::BlockDiagCholeskyPrecon<Type>::decoupledPrecondition
+(
+    Field<Type>& x,
+    const Field<Type>& b
+) const
+{
+    typedef CoeffField<Type> TypeCoeffField;
+
+    // Note: Assuming lower and upper triangle have the same active type
+
+    if (this->matrix_.symmetric())
+    {
+        const TypeCoeffField& UpperCoeff = this->matrix_.upper();
+
+        if (preconDiag_.activeType() == blockCoeffBase::SCALAR)
+        {
+            if (UpperCoeff.activeType() == blockCoeffBase::SCALAR)
+            {
+                ILUmultiply
+                (
+                    x,
+                    preconDiag_.asScalar(),
+                    UpperCoeff.asScalar(),
+                    b
+                );
+            }
+            else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR)
+            {
+                ILUmultiply
+                (
+                    x,
+                    preconDiag_.asScalar(),
+                    UpperCoeff.asLinear(),
+                    b
+                );
+            }
+        }
+        else if (preconDiag_.activeType() == blockCoeffBase::LINEAR)
+        {
+            if (UpperCoeff.activeType() == blockCoeffBase::SCALAR)
+            {
+                ILUmultiply
+                (
+                    x,
+                    preconDiag_.asLinear(),
+                    UpperCoeff.asScalar(),
+                    b
+                );
+            }
+            else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR)
+            {
+                ILUmultiply
+                (
+                    x,
+                    preconDiag_.asLinear(),
+                    UpperCoeff.asLinear(),
+                    b
+                );
+            }
+        }
+    }
+    else // Asymmetric matrix
+    {
+        const TypeCoeffField& LowerCoeff = this->matrix_.lower();
+        const TypeCoeffField& UpperCoeff = this->matrix_.upper();
+
+        if (preconDiag_.activeType() == blockCoeffBase::SCALAR)
+        {
+            if (UpperCoeff.activeType() == blockCoeffBase::SCALAR)
+            {
+                ILUmultiply
+                (
+                    x,
+                    preconDiag_.asScalar(),
+                    LowerCoeff.asScalar(),
+                    UpperCoeff.asScalar(),
+                    b
+                );
+            }
+            else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR)
+            {
+                ILUmultiply
+                (
+                    x,
+                    preconDiag_.asScalar(),
+                    LowerCoeff.asLinear(),
+                    UpperCoeff.asLinear(),
+                    b
+                );
+            }
+        }
+        else if (preconDiag_.activeType() == blockCoeffBase::LINEAR)
+        {
+            if (UpperCoeff.activeType() == blockCoeffBase::SCALAR)
+            {
+                ILUmultiply
+                (
+                    x,
+                    preconDiag_.asLinear(),
+                    LowerCoeff.asScalar(),
+                    UpperCoeff.asScalar(),
+                    b
+                );
+            }
+            else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR)
+            {
+                ILUmultiply
+                (
+                    x,
+                    preconDiag_.asLinear(),
+                    LowerCoeff.asLinear(),
+                    UpperCoeff.asLinear(),
+                    b
+                );
+            }
+        }
+    }
+}
+
+
+template<class Type>
+void Foam::BlockDiagCholeskyPrecon<Type>::decoupledPreconditionT
+(
+    Field<Type>& xT,
+    const Field<Type>& bT
+) const
+{
+    typedef CoeffField<Type> TypeCoeffField;
+
+    // Note: Assuming lower and upper triangle have the same active type
+
+    if (this->matrix_.symmetric())
+    {
+        precondition(xT, bT);
+    }
+    else // Asymmetric matrix
+    {
+        const TypeCoeffField& LowerCoeff = this->matrix_.lower();
+        const TypeCoeffField& UpperCoeff = this->matrix_.upper();
+
+        if (preconDiag_.activeType() == blockCoeffBase::SCALAR)
+        {
+            if (UpperCoeff.activeType() == blockCoeffBase::SCALAR)
+            {
+                ILUmultiplyTranspose
+                (
+                    xT,
+                    preconDiag_.asScalar(),
+                    LowerCoeff.asScalar(),
+                    UpperCoeff.asScalar(),
+                    bT
+                );
+            }
+            else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR)
+            {
+                ILUmultiplyTranspose
+                (
+                    xT,
+                    preconDiag_.asScalar(),
+                    LowerCoeff.asLinear(),
+                    UpperCoeff.asLinear(),
+                    bT
+                );
+            }
+        }
+        else if (preconDiag_.activeType() == blockCoeffBase::LINEAR)
+        {
+            if (UpperCoeff.activeType() == blockCoeffBase::SCALAR)
+            {
+                ILUmultiplyTranspose
+                (
+                    xT,
+                    preconDiag_.asLinear(),
+                    LowerCoeff.asScalar(),
+                    UpperCoeff.asScalar(),
+                    bT
+                );
+            }
+            else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR)
+            {
+                ILUmultiplyTranspose
+                (
+                    xT,
+                    preconDiag_.asLinear(),
+                    LowerCoeff.asLinear(),
+                    UpperCoeff.asLinear(),
+                    bT
+                );
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/blockDiagCholeskyPrecons.C

@@ -0,0 +1,42 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     3.2
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "blockLduMatrices.H"
+#include "blockLduPrecons.H"
+#include "blockDiagCholeskyPrecons.H"
+#include "addToRunTimeSelectionTable.H"
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+makeBlockPrecons(blockDiagCholeskyPrecon);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/blockDiagCholeskyPrecons.H

@@ -0,0 +1,66 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     3.2
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    BlockDiagCholeskyPrecon
+
+Description
+    Typedefs for Incomplete Cholesky preconditioning with no fill-in,
+    using the diagonal-of-diagonal for ILU decomposition.
+    Currently broken
+
+
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved.
+
+SourceFiles
+    blockDiagCholeskyPrecons.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef blockDiagCholeskyPrecons_H
+#define blockDiagCholeskyPrecons_H
+
+#include "scalarBlockDiagCholeskyPrecon.H"
+#include "tensorBlockDiagCholeskyPrecon.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+typedef BlockDiagCholeskyPrecon<scalar> blockDiagCholeskyPreconScalar;
+typedef BlockDiagCholeskyPrecon<vector> blockDiagCholeskyPreconVector;
+typedef BlockDiagCholeskyPrecon<tensor> blockDiagCholeskyPreconTensor;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/scalarBlockDiagCholeskyPrecon.C

@@ -0,0 +1,199 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     3.2
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef scalarBlockDiagCholeskyPrecon_H
+#define scalarBlockDiagCholeskyPrecon_H
+
+#include "BlockDiagCholeskyPrecon.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<>
+void Foam::BlockDiagCholeskyPrecon<scalar>::calcPreconDiag()
+{
+    // Precondition the diagonal
+    if (matrix_.symmetric())
+    {
+        const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr();
+        const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr();
+
+        // Get off-diagonal matrix coefficients
+        const scalarField& upper = matrix_.upper();
+
+        forAll (upper, coeffI)
+        {
+            preconDiag_[upperAddr[coeffI]] -=
+                sqr(upper[coeffI])/preconDiag_[lowerAddr[coeffI]];
+        }
+    }
+    else if (matrix_.asymmetric())
+    {
+        const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr();
+        const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr();
+
+        // Get off-diagonal matrix coefficients
+        const scalarField& upper = matrix_.upper();
+        const scalarField& lower = matrix_.lower();
+
+        forAll (upper, coeffI)
+        {
+            preconDiag_[upperAddr[coeffI]] -=
+                upper[coeffI]*lower[coeffI]/preconDiag_[lowerAddr[coeffI]];
+        }
+    }
+
+    // Invert the diagonal for future use
+    forAll (preconDiag_, i)
+    {
+        preconDiag_[i] = 1.0/preconDiag_[i];
+    }
+}
+
+
+template<>
+void Foam::BlockDiagCholeskyPrecon<scalar>::precondition
+(
+    scalarField& x,
+    const scalarField& b
+) const
+{
+    forAll(x, i)
+    {
+        x[i] = b[i]*preconDiag_[i];
+    }
+
+    if (matrix_.symmetric())
+    {
+        const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr();
+        const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr();
+
+        // Get off-diagonal matrix coefficients
+        const scalarField& upper = matrix_.upper();
+
+        forAll (upper, coeffI)
+        {
+            x[upperAddr[coeffI]] -=
+                preconDiag_[upperAddr[coeffI]]*
+                upper[coeffI]*x[lowerAddr[coeffI]];
+        }
+
+        forAllReverse (upper, coeffI)
+        {
+            x[lowerAddr[coeffI]] -=
+                preconDiag_[lowerAddr[coeffI]]*
+                upper[coeffI]*x[upperAddr[coeffI]];
+        }
+    }
+    else if (matrix_.asymmetric())
+    {
+        const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr();
+        const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr();
+        const unallocLabelList& losortAddr = matrix_.lduAddr().losortAddr();
+
+        // Get off-diagonal matrix coefficients
+        const scalarField& upper = matrix_.upper();
+        const scalarField& lower = matrix_.lower();
+
+        label losortCoeff;
+
+        forAll (lower, coeffI)
+        {
+            losortCoeff = losortAddr[coeffI];
+
+            x[upperAddr[losortCoeff]] -=
+                preconDiag_[upperAddr[losortCoeff]]*
+                lower[losortCoeff]*x[lowerAddr[losortCoeff]];
+        }
+
+        forAllReverse (upper, coeffI)
+        {
+            x[lowerAddr[coeffI]] -=
+                preconDiag_[lowerAddr[coeffI]]*
+                upper[coeffI]*x[upperAddr[coeffI]];
+        }
+    }
+}
+
+
+template<>
+void Foam::BlockDiagCholeskyPrecon<scalar>::preconditionT
+(
+    scalarField& xT,
+    const scalarField& bT
+) const
+{
+    if (matrix_.symmetric())
+    {
+        precondition(xT, bT);
+    }
+
+    forAll(xT, i)
+    {
+        xT[i] = bT[i]*preconDiag_[i];
+    }
+
+    if (matrix_.asymmetric())
+    {
+        const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr();
+        const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr();
+        const unallocLabelList& losortAddr = matrix_.lduAddr().losortAddr();
+
+        // Get off-diagonal matrix coefficients
+        const scalarField& upper = matrix_.upper();
+        const scalarField& lower = matrix_.lower();
+
+        label losortCoeff;
+
+        forAll (lower, coeffI)
+        {
+            xT[upperAddr[coeffI]] -=
+                preconDiag_[upperAddr[coeffI]]*
+                upper[coeffI]*xT[lowerAddr[coeffI]];
+        }
+
+        forAllReverse (upper, coeffI)
+        {
+            losortCoeff = losortAddr[coeffI];
+
+            xT[lowerAddr[losortCoeff]] -=
+                preconDiag_[lowerAddr[losortCoeff]]*
+                lower[losortCoeff]*xT[upperAddr[losortCoeff]];
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/scalarBlockDiagCholeskyPrecon.H

@@ -0,0 +1,75 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     3.2
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    BlockDiagCholeskyPrecon
+
+Description
+    Template specialisation for scalar block Cholesky preconditioning
+
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved.
+
+SourceFiles
+    scalarBlockDiagCholeskyPrecon.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef scalarBlockDiagCholeskyPrecon_H
+#define scalarBlockDiagCholeskyPrecon_H
+
+#include "BlockDiagCholeskyPrecon.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<>
+void Foam::BlockDiagCholeskyPrecon<scalar>::calcPreconDiag();
+
+template<>
+void Foam::BlockDiagCholeskyPrecon<scalar>::precondition
+(
+    scalarField& x,
+    const scalarField& b
+) const;
+
+
+template<>
+void Foam::BlockDiagCholeskyPrecon<scalar>::preconditionT
+(
+    scalarField& xT,
+    const scalarField& bT
+) const;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/tensorBlockDiagCholeskyPrecon.C

@@ -0,0 +1,76 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     3.2
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef tensorBlockDiagCholeskyPrecon_H
+#define tensorBlockDiagCholeskyPrecon_H
+
+#include "BlockDiagCholeskyPrecon.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<>
+void Foam::BlockDiagCholeskyPrecon<tensor>::calcPreconDiag()
+{
+    // Decoupled version
+    calcDecoupledPreconDiag();
+}
+
+
+template<>
+void Foam::BlockDiagCholeskyPrecon<tensor>::precondition
+(
+    tensorField& x,
+    const tensorField& b
+) const
+{
+    // Decoupled version
+    decoupledPrecondition(x, b);
+}
+
+
+template<>
+void Foam::BlockDiagCholeskyPrecon<tensor>::preconditionT
+(
+    tensorField& xT,
+    const tensorField& bT
+) const
+{
+    // Decoupled version
+    decoupledPreconditionT(xT, bT);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/tensorBlockDiagCholeskyPrecon.H

@@ -0,0 +1,75 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     3.2
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    BlockDiagCholeskyPrecon
+
+Description
+    Template specialisation for tensor block Cholesky preconditioning
+
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved.
+
+SourceFiles
+    tensorBlockDiagCholeskyPrecon.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef tensorBlockDiagCholeskyPrecon_H
+#define tensorBlockDiagCholeskyPrecon_H
+
+#include "BlockDiagCholeskyPrecon.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<>
+void Foam::BlockDiagCholeskyPrecon<tensor>::calcPreconDiag();
+
+template<>
+void Foam::BlockDiagCholeskyPrecon<tensor>::precondition
+(
+    tensorField& x,
+    const tensorField& b
+) const;
+
+
+template<>
+void Foam::BlockDiagCholeskyPrecon<tensor>::preconditionT
+(
+    tensorField& xT,
+    const tensorField& bT
+) const;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.C

@@ -34,6 +34,79 @@ Author
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+template<class Type>
+void Foam::BlockGaussSeidelPrecon<Type>::calcInvDiag()
+{
+    typedef CoeffField<Type> TypeCoeffField;
+    typedef typename TypeCoeffField::linearTypeField linearTypeField;
+    typedef typename TypeCoeffField::linearType valueType;
+
+    typedef typename TypeCoeffField::squareTypeField squareTypeField;
+
+    // Note: Cannot use inv since the square coefficient requires
+    // special treatment.  HJ, 20/Aug/2015
+
+    // Get reference to diagonal
+    const TypeCoeffField& d = this->matrix_.diag();
+
+    if (d.activeType() == blockCoeffBase::SCALAR)
+    {
+        invDiag_.asScalar() = 1/d.asScalar();
+    }
+    else if (d.activeType() == blockCoeffBase::LINEAR)
+    {
+        invDiag_.asLinear() =
+            cmptDivide
+            (
+                linearTypeField(d.size(), pTraits<valueType>::one),
+                d.asLinear()
+            );
+    }
+    else if (d.activeType() == blockCoeffBase::SQUARE)
+    {
+        // For square diagonal invert diagonal only and store the rest
+        // info LUDiag coefficient.  This avoids full inverse of the
+        // diagonal matrix.  HJ, 20/Aug/2015
+        Info<< "Square diag inverse" << endl;
+
+        // Get reference to active diag
+        const squareTypeField& activeDiag = d.asSquare();
+
+        // Get reference to LU: remove diagonal from active diag
+        squareTypeField& luDiag = LUDiag_.asSquare();
+
+        linearTypeField lf(activeDiag.size());
+
+        // Take out the diagonal from the diag as a linear type
+        contractLinear(lf, activeDiag);
+
+        // Expand diagonal only to full square type and store into luDiag
+        expandLinear(luDiag, lf);
+
+        // Keep only off-diagonal in ldDiag.
+        // Note change of sign to avoid multiplication with -1 when moving
+        // to the other side.  HJ, 20/Aug/2015
+        luDiag -= activeDiag;
+
+        // Inverse is the inverse of diagonal only
+        invDiag_.asLinear() =
+            cmptDivide
+            (
+                linearTypeField(lf.size(), pTraits<valueType>::one),
+                lf
+            );
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "void BlockGaussSeidelPrecon<Type>::calcInvDiag()"
+        )   << "Problem with coefficient type morphing."
+            << abort(FatalError);
+    }
+}
+
+
 // Block sweep, symmetric matrix
 template<class Type>
 template<class DiagType, class ULType>
@@ -46,7 +119,8 @@ void Foam::BlockGaussSeidelPrecon<Type>::BlockSweep
 ) const
 {
     const unallocLabelList& u = this->matrix_.lduAddr().upperAddr();
-    const unallocLabelList& ownStart = this->matrix_.lduAddr().ownerStartAddr();
+    const unallocLabelList& ownStart =
+        this->matrix_.lduAddr().ownerStartAddr();
 
     const label nRows = ownStart.size() - 1;
 
@@ -155,7 +229,8 @@ void Foam::BlockGaussSeidelPrecon<Type>::BlockSweep
 ) const
 {
     const unallocLabelList& u = this->matrix_.lduAddr().upperAddr();
-    const unallocLabelList& ownStart = this->matrix_.lduAddr().ownerStartAddr();
+    const unallocLabelList& ownStart =
+        this->matrix_.lduAddr().ownerStartAddr();
 
     const label nRows = ownStart.size() - 1;
 
@@ -247,6 +322,43 @@ void Foam::BlockGaussSeidelPrecon<Type>::BlockSweep
 }
 
 
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::BlockGaussSeidelPrecon<Type>::BlockGaussSeidelPrecon
+(
+    const BlockLduMatrix<Type>& matrix
+)
+:
+    BlockLduPrecon<Type>(matrix),
+    invDiag_(matrix.lduAddr().size()),
+    LUDiag_(matrix.lduAddr().size()),
+    bPlusLU_(),
+    bPrime_(matrix.lduAddr().size()),
+    nSweeps_(1)
+{
+    calcInvDiag();
+}
+
+
+template<class Type>
+Foam::BlockGaussSeidelPrecon<Type>::BlockGaussSeidelPrecon
+(
+    const BlockLduMatrix<Type>& matrix,
+    const dictionary& dict
+)
+:
+    BlockLduPrecon<Type>(matrix),
+    invDiag_(matrix.lduAddr().size()),
+    LUDiag_(matrix.lduAddr().size()),
+    bPlusLU_(),
+    bPrime_(matrix.lduAddr().size()),
+    nSweeps_(readLabel(dict.lookup("nSweeps")))
+{
+    calcInvDiag();
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
@@ -260,14 +372,11 @@ void Foam::BlockGaussSeidelPrecon<Type>::precondition
 
     if (this->matrix_.diagonal())
     {
-        TypeCoeffField dDCoeff = inv(this->matrix_.diag());
-
-        multiply(x, dDCoeff, b);
+        multiply(x, invDiag_, b);
     }
     else if (this->matrix_.symmetric())
     {
-        TypeCoeffField dDCoeff = inv(this->matrix_.diag());
-
+        const TypeCoeffField& DiagCoeff = this->matrix_.diag();
         const TypeCoeffField& UpperCoeff = this->matrix_.upper();
 
         // Note
@@ -281,14 +390,14 @@ void Foam::BlockGaussSeidelPrecon<Type>::precondition
 
         // Note: Assuming lower and upper triangle have the same active type
 
-        if (dDCoeff.activeType() == blockCoeffBase::SCALAR)
+        if (DiagCoeff.activeType() == blockCoeffBase::SCALAR)
         {
             if (UpperCoeff.activeType() == blockCoeffBase::SCALAR)
             {
                 BlockSweep
                 (
                     x,
-                    dDCoeff.asScalar(),
+                    invDiag_.asScalar(),
                     UpperCoeff.asScalar(),
                     b
                 );
@@ -298,7 +407,7 @@ void Foam::BlockGaussSeidelPrecon<Type>::precondition
                 BlockSweep
                 (
                     x,
-                    dDCoeff.asScalar(),
+                    invDiag_.asScalar(),
                     UpperCoeff.asLinear(),
                     b
                 );
@@ -308,20 +417,20 @@ void Foam::BlockGaussSeidelPrecon<Type>::precondition
                 BlockSweep
                 (
                     x,
-                    dDCoeff.asScalar(),
+                    invDiag_.asScalar(),
                     UpperCoeff.asSquare(),
                     b
                 );
             }
         }
-        else if (dDCoeff.activeType() == blockCoeffBase::LINEAR)
+        else if (DiagCoeff.activeType() == blockCoeffBase::LINEAR)
         {
             if (UpperCoeff.activeType() == blockCoeffBase::SCALAR)
             {
                 BlockSweep
                 (
                     x,
-                    dDCoeff.asLinear(),
+                    invDiag_.asLinear(),
                     UpperCoeff.asScalar(),
                     b
                 );
@@ -331,7 +440,7 @@ void Foam::BlockGaussSeidelPrecon<Type>::precondition
                 BlockSweep
                 (
                     x,
-                    dDCoeff.asLinear(),
+                    invDiag_.asLinear(),
                     UpperCoeff.asLinear(),
                     b
                 );
@@ -341,42 +450,57 @@ void Foam::BlockGaussSeidelPrecon<Type>::precondition
                 BlockSweep
                 (
                     x,
-                    dDCoeff.asLinear(),
+                    invDiag_.asLinear(),
                     UpperCoeff.asSquare(),
                     b
                 );
             }
         }
-        else if (dDCoeff.activeType() == blockCoeffBase::SQUARE)
+        else if (DiagCoeff.activeType() == blockCoeffBase::SQUARE)
         {
+            // Add diag coupling to b
+            if (bPlusLU_.empty())
+            {
+                bPlusLU_.setSize(b.size(), pTraits<Type>::zero);
+            }
+
+            // Multiply overwrites bPlusLU_: no need to initialise
+            // Change of sign accounted via change of sign of bPlusLU_
+            // HJ, 20/Aug/2015
+            multiply(bPlusLU_, LUDiag_, x);
+            bPlusLU_ += b;
+
             if (UpperCoeff.activeType() == blockCoeffBase::SCALAR)
             {
+                // Note linear diag inversed due to decoupling
                 BlockSweep
                 (
                     x,
-                    dDCoeff.asSquare(),
+                    invDiag_.asLinear(),
                     UpperCoeff.asScalar(),
-                    b
+                    bPlusLU_
                 );
             }
             else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR)
             {
+                // Note linear diag inversed due to decoupling
                 BlockSweep
                 (
                     x,
-                    dDCoeff.asSquare(),
+                    invDiag_.asLinear(),
                     UpperCoeff.asLinear(),
-                    b
+                    bPlusLU_
                 );
             }
             else if (UpperCoeff.activeType() == blockCoeffBase::SQUARE)
             {
+                // Note linear diag inversed due to decoupling
                 BlockSweep
                 (
                     x,
-                    dDCoeff.asSquare(),
+                    invDiag_.asLinear(),
                     UpperCoeff.asSquare(),
-                    b
+                    bPlusLU_
                 );
             }
         }
@@ -395,8 +519,7 @@ void Foam::BlockGaussSeidelPrecon<Type>::precondition
     }
     else if (this->matrix_.asymmetric())
     {
-        TypeCoeffField dDCoeff = inv(this->matrix_.diag());
-
+        const TypeCoeffField& DiagCoeff = this->matrix_.diag();
         const TypeCoeffField& LowerCoeff = this->matrix_.lower();
         const TypeCoeffField& UpperCoeff = this->matrix_.upper();
 
@@ -411,14 +534,14 @@ void Foam::BlockGaussSeidelPrecon<Type>::precondition
 
         //Note: Assuming lower and upper triangle have the same active type
 
-        if (dDCoeff.activeType() == blockCoeffBase::SCALAR)
+        if (DiagCoeff.activeType() == blockCoeffBase::SCALAR)
         {
             if (UpperCoeff.activeType() == blockCoeffBase::SCALAR)
             {
                 BlockSweep
                 (
                     x,
-                    dDCoeff.asScalar(),
+                    invDiag_.asScalar(),
                     LowerCoeff.asScalar(),
                     UpperCoeff.asScalar(),
                     b
@@ -429,7 +552,7 @@ void Foam::BlockGaussSeidelPrecon<Type>::precondition
                 BlockSweep
                 (
                     x,
-                    dDCoeff.asScalar(),
+                    invDiag_.asScalar(),
                     LowerCoeff.asLinear(),
                     UpperCoeff.asLinear(),
                     b
@@ -440,21 +563,21 @@ void Foam::BlockGaussSeidelPrecon<Type>::precondition
                 BlockSweep
                 (
                     x,
-                    dDCoeff.asScalar(),
+                    invDiag_.asScalar(),
                     LowerCoeff.asSquare(),
                     UpperCoeff.asSquare(),
                     b
                 );
             }
         }
-        else if (dDCoeff.activeType() == blockCoeffBase::LINEAR)
+        else if (DiagCoeff.activeType() == blockCoeffBase::LINEAR)
         {
             if (UpperCoeff.activeType() == blockCoeffBase::SCALAR)
             {
                 BlockSweep
                 (
                     x,
-                    dDCoeff.asLinear(),
+                    invDiag_.asLinear(),
                     LowerCoeff.asScalar(),
                     UpperCoeff.asScalar(),
                     b
@@ -465,7 +588,7 @@ void Foam::BlockGaussSeidelPrecon<Type>::precondition
                 BlockSweep
                 (
                     x,
-                    dDCoeff.asLinear(),
+                    invDiag_.asLinear(),
                     LowerCoeff.asLinear(),
                     UpperCoeff.asLinear(),
                     b
@@ -476,46 +599,61 @@ void Foam::BlockGaussSeidelPrecon<Type>::precondition
                 BlockSweep
                 (
                     x,
-                    dDCoeff.asLinear(),
+                    invDiag_.asLinear(),
                     LowerCoeff.asSquare(),
                     UpperCoeff.asSquare(),
                     b
                 );
             }
         }
-        else if (dDCoeff.activeType() == blockCoeffBase::SQUARE)
+        else if (DiagCoeff.activeType() == blockCoeffBase::SQUARE)
         {
+            // Add diag coupling to b
+            if (bPlusLU_.empty())
+            {
+                bPlusLU_.setSize(b.size(), pTraits<Type>::zero);
+            }
+
+            // Multiply overwrites bPlusLU_: no need to initialise
+            // Change of sign accounted via change of sign of bPlusLU_
+            // HJ, 20/Aug/2015
+            multiply(bPlusLU_, LUDiag_, x);
+            bPlusLU_ += b;
+
             if (UpperCoeff.activeType() == blockCoeffBase::SCALAR)
             {
+                // Note linear diag inversed due to decoupling
                 BlockSweep
                 (
                     x,
-                    dDCoeff.asSquare(),
+                    invDiag_.asLinear(),
                     LowerCoeff.asScalar(),
                     UpperCoeff.asScalar(),
-                    b
+                    bPlusLU_
                 );
             }
             else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR)
             {
+                // Note linear diag inversed due to decoupling
                 BlockSweep
                 (
                     x,
-                    dDCoeff.asSquare(),
+                    invDiag_.asLinear(),
                     LowerCoeff.asLinear(),
                     UpperCoeff.asLinear(),
-                    b
+                    bPlusLU_
                 );
             }
             else if (UpperCoeff.activeType() == blockCoeffBase::SQUARE)
             {
+                // Note linear diag inversed due to decoupling
                 BlockSweep
                 (
                     x,
-                    dDCoeff.asSquare(),
+                    invDiag_.asLinear(),
                     LowerCoeff.asSquare(),
                     UpperCoeff.asSquare(),
-                    b
+                    bPlusLU_
                 );
             }
         }
@@ -558,14 +696,11 @@ void Foam::BlockGaussSeidelPrecon<Type>::preconditionT
 
     if (this->matrix_.diagonal())
     {
-        TypeCoeffField dDCoeff = inv(this->matrix_.diag());
-
-        multiply(xT, dDCoeff, bT);
+        multiply(xT, invDiag_, bT);
     }
     else if (this->matrix_.symmetric() || this->matrix_.asymmetric())
     {
-        TypeCoeffField dDCoeff = inv(this->matrix_.diag());
-
+        const TypeCoeffField& DiagCoeff = this->matrix_.lower();
         const TypeCoeffField& LowerCoeff = this->matrix_.lower();
         const TypeCoeffField& UpperCoeff = this->matrix_.upper();
 
@@ -580,7 +715,7 @@ void Foam::BlockGaussSeidelPrecon<Type>::preconditionT
 
         //Note: Assuming lower and upper triangle have the same active type
 
-        if (dDCoeff.activeType() == blockCoeffBase::SCALAR)
+        if (DiagCoeff.activeType() == blockCoeffBase::SCALAR)
         {
             if (UpperCoeff.activeType() == blockCoeffBase::SCALAR)
             {
@@ -588,7 +723,7 @@ void Foam::BlockGaussSeidelPrecon<Type>::preconditionT
                 BlockSweep
                 (
                     xT,
-                    dDCoeff.asScalar(),
+                    invDiag_.asScalar(),
                     UpperCoeff.asScalar(),
                     LowerCoeff.asScalar(),
                     bT
@@ -600,7 +735,7 @@ void Foam::BlockGaussSeidelPrecon<Type>::preconditionT
                 BlockSweep
                 (
                     xT,
-                    dDCoeff.asScalar(),
+                    invDiag_.asScalar(),
                     UpperCoeff.asLinear(),
                     LowerCoeff.asLinear(),
                     bT
@@ -612,14 +747,14 @@ void Foam::BlockGaussSeidelPrecon<Type>::preconditionT
                 BlockSweep
                 (
                     xT,
-                    dDCoeff.asScalar(),
+                    invDiag_.asScalar(),
                     UpperCoeff.asSquare(),
                     LowerCoeff.asSquare(),
                     bT
                 );
             }
         }
-        else if (dDCoeff.activeType() == blockCoeffBase::LINEAR)
+        else if (DiagCoeff.activeType() == blockCoeffBase::LINEAR)
         {
             if (UpperCoeff.activeType() == blockCoeffBase::SCALAR)
             {
@@ -627,7 +762,7 @@ void Foam::BlockGaussSeidelPrecon<Type>::preconditionT
                 BlockSweep
                 (
                     xT,
-                    dDCoeff.asLinear(),
+                    invDiag_.asLinear(),
                     UpperCoeff.asScalar(),
                     LowerCoeff.asScalar(),
                     bT
@@ -639,7 +774,7 @@ void Foam::BlockGaussSeidelPrecon<Type>::preconditionT
                 BlockSweep
                 (
                     xT,
-                    dDCoeff.asLinear(),
+                    invDiag_.asLinear(),
                     UpperCoeff.asLinear(),
                     LowerCoeff.asLinear(),
                     bT
@@ -651,49 +786,64 @@ void Foam::BlockGaussSeidelPrecon<Type>::preconditionT
                 BlockSweep
                 (
                     xT,
-                    dDCoeff.asLinear(),
+                    invDiag_.asLinear(),
                     UpperCoeff.asSquare(),
                     LowerCoeff.asSquare(),
                     bT
                 );
             }
         }
-        else if (dDCoeff.activeType() == blockCoeffBase::SQUARE)
+        else if (DiagCoeff.activeType() == blockCoeffBase::SQUARE)
         {
+            // Add diag coupling to b
+            if (bPlusLU_.empty())
+            {
+                bPlusLU_.setSize(bT.size(), pTraits<Type>::zero);
+            }
+
+            // Multiply overwrites bPlusLU_: no need to initialise
+            // Change of sign accounted via change of sign of bPlusLU_
+            // HJ, 20/Aug/2015
+            multiply(bPlusLU_, LUDiag_, xT);
+            bPlusLU_ += bT;
+
             if (UpperCoeff.activeType() == blockCoeffBase::SCALAR)
             {
                 // Transpose multiplication - swap lower and upper coeff arrays
+                // Note linear diag inversed due to decoupling
                 BlockSweep
                 (
                     xT,
-                    dDCoeff.asSquare(),
+                    invDiag_.asLinear(),
                     UpperCoeff.asScalar(),
                     LowerCoeff.asScalar(),
-                    bT
+                    bPlusLU_
                 );
             }
             else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR)
             {
                 // Transpose multiplication - swap lower and upper coeff arrays
+                // Note linear diag inversed due to decoupling
                 BlockSweep
                 (
                     xT,
-                    dDCoeff.asSquare(),
+                    invDiag_.asLinear(),
                     UpperCoeff.asLinear(),
                     LowerCoeff.asLinear(),
-                    bT
+                    bPlusLU_
                 );
             }
             else if (UpperCoeff.activeType() == blockCoeffBase::SQUARE)
             {
                 // Transpose multiplication - swap lower and upper coeff arrays
+                // Note linear diag inversed due to decoupling
                 BlockSweep
                 (
                     xT,
-                    dDCoeff.asSquare(),
+                    invDiag_.asLinear(),
                     UpperCoeff.asSquare(),
                     LowerCoeff.asSquare(),
-                    bT
+                    bPlusLU_
                 );
             }
         }

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.H

@@ -58,6 +58,16 @@ class BlockGaussSeidelPrecon
 {
     // Private Data
 
+        //- Inverse of diagonal diagonal
+        CoeffField<Type> invDiag_;
+
+        //- Off-diag part of diagonal
+        CoeffField<Type> LUDiag_;
+
+        //- Temporary space for updated decoupled source
+        //  Initialised with zero size and resized on first use
+        mutable Field<Type> bPlusLU_;
+
         //- Temporary space for solution intermediate
         mutable Field<Type> bPrime_;
 
@@ -74,6 +84,9 @@ class BlockGaussSeidelPrecon
         void operator=(const BlockGaussSeidelPrecon&);
 
 
+        //- Calculate inverse diagonal
+        void calcInvDiag();
+
         // Block Gauss-Seidel sweep, symetric matrix
         template<class DiagType, class ULType>
         void BlockSweep
@@ -98,6 +111,9 @@ class BlockGaussSeidelPrecon
 
         // Decoupled operations, used in template specialisation
 
+           //- Calculate inverse diagonal, decoupled version
+           void calcDecoupledInvDiag();
+
             //- Execute preconditioning, decoupled version
             void decoupledPrecondition
             (
@@ -126,24 +142,14 @@ public:
         BlockGaussSeidelPrecon
         (
             const BlockLduMatrix<Type>& matrix
-        )
-        :
-            BlockLduPrecon<Type>(matrix),
-            bPrime_(matrix.lduAddr().size()),
-            nSweeps_(1)
-        {}
+        );
 
         //- Construct from components
         BlockGaussSeidelPrecon
         (
             const BlockLduMatrix<Type>& matrix,
             const dictionary& dict
-        )
-        :
-            BlockLduPrecon<Type>(matrix),
-            bPrime_(matrix.lduAddr().size()),
-            nSweeps_(readLabel(dict.lookup("nSweeps")))
-        {}
+        );
 
 
     // Destructor

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPreconDecoupled.C

@@ -34,6 +34,19 @@ Author
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+template<class Type>
+void Foam::BlockGaussSeidelPrecon<Type>::calcDecoupledInvDiag()
+{
+    // Get reference to diagonal and obtain inverse by casting
+    typedef CoeffField<Type> TypeCoeffField;
+
+    const TypeCoeffField& d = this->matrix_.diag();
+    const DecoupledCoeffField<Type>& dd = d;
+
+    invDiag_ = CoeffField<Type>(inv(dd)());
+}
+
+
 template<class Type>
 void Foam::BlockGaussSeidelPrecon<Type>::decoupledPrecondition
 (
@@ -45,14 +58,17 @@ void Foam::BlockGaussSeidelPrecon<Type>::decoupledPrecondition
 
     if (this->matrix_.diagonal())
     {
-        TypeCoeffField dDCoeff = inv(this->matrix_.diag());
-
-        multiply(x, dDCoeff, b);
+        if (invDiag_.activeType() == blockCoeffBase::SCALAR)
+        {
+            x = invDiag_.asScalar()*b;
+        }
+        else if (invDiag_.activeType() == blockCoeffBase::LINEAR)
+        {
+            x = cmptMultiply(invDiag_.asLinear(), b);
+        }
     }
     else if (this->matrix_.symmetric() || this->matrix_.asymmetric())
     {
-        TypeCoeffField dDCoeff = inv(this->matrix_.diag());
-
         const TypeCoeffField& LowerCoeff = this->matrix_.lower();
         const TypeCoeffField& UpperCoeff = this->matrix_.upper();
 
@@ -67,14 +83,14 @@ void Foam::BlockGaussSeidelPrecon<Type>::decoupledPrecondition
 
         // Note: Assuming lower and upper triangle have the same active type
 
-        if (dDCoeff.activeType() == blockCoeffBase::SCALAR)
+        if (invDiag_.activeType() == blockCoeffBase::SCALAR)
         {
             if (UpperCoeff.activeType() == blockCoeffBase::SCALAR)
             {
                 BlockSweep
                 (
                     x,
-                    dDCoeff.asScalar(),
+                    invDiag_.asScalar(),
                     LowerCoeff.asScalar(),
                     UpperCoeff.asScalar(),
                     b
@@ -85,21 +101,21 @@ void Foam::BlockGaussSeidelPrecon<Type>::decoupledPrecondition
                 BlockSweep
                 (
                     x,
-                    dDCoeff.asScalar(),
+                    invDiag_.asScalar(),
                     LowerCoeff.asLinear(),
                     UpperCoeff.asLinear(),
                     b
                 );
             }
         }
-        else if (dDCoeff.activeType() == blockCoeffBase::LINEAR)
+        else if (invDiag_.activeType() == blockCoeffBase::LINEAR)
         {
             if (UpperCoeff.activeType() == blockCoeffBase::SCALAR)
             {
                 BlockSweep
                 (
                     x,
-                    dDCoeff.asLinear(),
+                    invDiag_.asLinear(),
                     LowerCoeff.asScalar(),
                     UpperCoeff.asScalar(),
                     b
@@ -110,7 +126,7 @@ void Foam::BlockGaussSeidelPrecon<Type>::decoupledPrecondition
                 BlockSweep
                 (
                     x,
-                    dDCoeff.asLinear(),
+                    invDiag_.asLinear(),
                     LowerCoeff.asLinear(),
                     UpperCoeff.asLinear(),
                     b
@@ -156,14 +172,17 @@ void Foam::BlockGaussSeidelPrecon<Type>::decoupledPreconditionT
 
     if (this->matrix_.diagonal())
     {
-        TypeCoeffField dDCoeff = inv(this->matrix_.diag());
-
-        multiply(xT, dDCoeff, bT);
+        if (invDiag_.activeType() == blockCoeffBase::SCALAR)
+        {
+            xT = invDiag_.asScalar()*bT;
+        }
+        else if (invDiag_.activeType() == blockCoeffBase::LINEAR)
+        {
+            xT = cmptMultiply(invDiag_.asLinear(), bT);
+        }
     }
     else if (this->matrix_.symmetric() || this->matrix_.asymmetric())
     {
-        TypeCoeffField dDCoeff = inv(this->matrix_.diag());
-
         const TypeCoeffField& LowerCoeff = this->matrix_.lower();
         const TypeCoeffField& UpperCoeff = this->matrix_.upper();
 
@@ -178,7 +197,7 @@ void Foam::BlockGaussSeidelPrecon<Type>::decoupledPreconditionT
 
         // Note: Assuming lower and upper triangle have the same active type
 
-        if (dDCoeff.activeType() == blockCoeffBase::SCALAR)
+        if (invDiag_.activeType() == blockCoeffBase::SCALAR)
         {
             if (UpperCoeff.activeType() == blockCoeffBase::SCALAR)
             {
@@ -186,7 +205,7 @@ void Foam::BlockGaussSeidelPrecon<Type>::decoupledPreconditionT
                 BlockSweep
                 (
                     xT,
-                    dDCoeff.asScalar(),
+                    invDiag_.asScalar(),
                     UpperCoeff.asScalar(),
                     LowerCoeff.asScalar(),
                     bT
@@ -198,14 +217,14 @@ void Foam::BlockGaussSeidelPrecon<Type>::decoupledPreconditionT
                 BlockSweep
                 (
                     xT,
-                    dDCoeff.asScalar(),
+                    invDiag_.asScalar(),
                     UpperCoeff.asLinear(),
                     LowerCoeff.asLinear(),
                     bT
                 );
             }
         }
-        else if (dDCoeff.activeType() == blockCoeffBase::LINEAR)
+        else if (invDiag_.activeType() == blockCoeffBase::LINEAR)
         {
             if (UpperCoeff.activeType() == blockCoeffBase::SCALAR)
             {
@@ -213,7 +232,7 @@ void Foam::BlockGaussSeidelPrecon<Type>::decoupledPreconditionT
                 BlockSweep
                 (
                     xT,
-                    dDCoeff.asLinear(),
+                    invDiag_.asLinear(),
                     UpperCoeff.asScalar(),
                     LowerCoeff.asScalar(),
                     bT
@@ -225,7 +244,7 @@ void Foam::BlockGaussSeidelPrecon<Type>::decoupledPreconditionT
                 BlockSweep
                 (
                     xT,
-                    dDCoeff.asLinear(),
+                    invDiag_.asLinear(),
                     UpperCoeff.asLinear(),
                     LowerCoeff.asLinear(),
                     bT

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/scalarBlockGaussSeidelPrecon.C

@@ -44,38 +44,11 @@ namespace Foam
 {
 
 template<>
-template<>
-void BlockGaussSeidelPrecon<scalar>::BlockSweep
-(
-    Field<scalar>& x,
-    const Field<scalar>& dD,
-    const Field<scalar>& upper,
-    const Field<scalar>& b
-) const
+void BlockGaussSeidelPrecon<scalar>::calcInvDiag()
 {
-    FatalErrorIn
-    (
-        "BlockGaussSeidelPrecon<scalar>::BlockSweep(...)"
-    )   << "Function not implemented for Type=scalar. " << endl
-        << abort(FatalError);
-}
-
-template<>
-template<>
-void BlockGaussSeidelPrecon<scalar>::BlockSweep
-(
-    Field<scalar>& x,
-    const Field<scalar>& dD,
-    const Field<scalar>& upper,
-    const Field<scalar>& lower,
-    const Field<scalar>& b
-) const
-{
-    FatalErrorIn
-    (
-        "BlockGaussSeidelPrecon<scalar>::BlockSweep(...)"
-    )   << "Function not implemented for Type=scalar. " << endl
-        << abort(FatalError);
+    // Direct inversion of diagonal is sufficient, as the diagonal
+    // is linear.  HJ, 20/Aug/2015
+    invDiag_ = inv(this->matrix_.diag());
 }
 
 

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/scalarBlockGaussSeidelPrecon.H

@@ -43,25 +43,7 @@ namespace Foam
 {
 
 template<>
-template<>
-void BlockGaussSeidelPrecon<scalar>::BlockSweep
-(
-    scalarField& x,
-    const scalarField& dD,
-    const scalarField& upper,
-    const scalarField& b
-) const;
-
-template<>
-template<>
-void BlockGaussSeidelPrecon<scalar>::BlockSweep
-(
-    scalarField& x,
-    const scalarField& dD,
-    const scalarField& upper,
-    const scalarField& lower,
-    const scalarField& b
-) const;
+void BlockGaussSeidelPrecon<scalar>::calcInvDiag();
 
 
 template<>

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/tensorBlockGaussSeidelPrecon.C

@@ -47,39 +47,9 @@ namespace Foam
 {
 
 template<>
-template<>
-void BlockGaussSeidelPrecon<tensor>::BlockSweep
-(
-    Field<tensor>& x,
-    const Field<tensor>& dD,
-    const Field<tensor>& upper,
-    const Field<tensor>& b
-) const
+void BlockGaussSeidelPrecon<tensor>::calcInvDiag()
 {
-    FatalErrorIn
-    (
-        "Foam::BlockGaussSeidelPrecon<tensor>::BlockSweep(...)"
-    )   << "Function not implemented for Type=tensor. " << endl
-        << abort(FatalError);
-}
-
-
-template<>
-template<>
-void BlockGaussSeidelPrecon<tensor>::BlockSweep
-(
-    Field<tensor>& x,
-    const Field<tensor>& dD,
-    const Field<tensor>& upper,
-    const Field<tensor>& lower,
-    const Field<tensor>& b
-) const
-{
-    FatalErrorIn
-    (
-        "Foam::BlockGaussSeidelPrecon<tensor>::BlockSweep(...)"
-    )   << "Function not implemented for Type=tensor. " << endl
-        << abort(FatalError);
+    calcDecoupledInvDiag();
 }
 
 
@@ -91,7 +61,7 @@ void BlockGaussSeidelPrecon<tensor>::precondition
 ) const
 {
     // Decoupled version
-//     decoupledPrecondition(x, b);
+    decoupledPrecondition(x, b);
 }
 
 
@@ -103,7 +73,7 @@ void BlockGaussSeidelPrecon<tensor>::preconditionT
 ) const
 {
     // Decoupled version
-//     decoupledPreconditionT(xT, bT);
+    decoupledPreconditionT(xT, bT);
 }
 
 

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/tensorBlockGaussSeidelPrecon.H

@@ -43,26 +43,7 @@ namespace Foam
 {
 
 template<>
-template<>
-void BlockGaussSeidelPrecon<tensor>::BlockSweep
-(
-    Field<tensor>& x,
-    const Field<tensor>& dD,
-    const Field<tensor>& upper,
-    const Field<tensor>& b
-) const;
-
-
-template<>
-template<>
-void BlockGaussSeidelPrecon<tensor>::BlockSweep
-(
-    Field<tensor>& x,
-    const Field<tensor>& dD,
-    const Field<tensor>& upper,
-    const Field<tensor>& lower,
-    const Field<tensor>& b
-) const;
+void BlockGaussSeidelPrecon<tensor>::calcInvDiag();
 
 
 template<>

src/foam/primitives/VectorN/DiagTensorNI.H

@@ -40,6 +40,7 @@ const DiagTensorN<Cmpt, length> DiagTensorN<Cmpt, length>::zero(0);
 template <class Cmpt, int length>
 const DiagTensorN<Cmpt, length> DiagTensorN<Cmpt, length>::one(1);
 
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct null
@@ -63,7 +64,12 @@ inline DiagTensorN<Cmpt, length>::DiagTensorN
 template <class Cmpt, int length>
 inline DiagTensorN<Cmpt, length>::DiagTensorN(const Cmpt& tx)
 {
-    VectorSpaceOps<DiagTensorN<Cmpt, length>::nComponents,0>::eqOpS(*this, tx, eqOp<Cmpt>());
+    VectorSpaceOps<DiagTensorN<Cmpt, length>::nComponents,0>::eqOpS
+    (
+        *this,
+        tx,
+        eqOp<Cmpt>()
+    );
 }
 
 
@@ -94,10 +100,18 @@ inline DiagTensorN<Cmpt, length> DiagTensorN<Cmpt, length>::T() const
 
 //- Assign to a SphericalTensorN
 template <class Cmpt, int length>
-inline void DiagTensorN<Cmpt, length>::operator=(const SphericalTensorN<Cmpt, length>& st)
+inline void DiagTensorN<Cmpt, length>::operator=
+(
+    const SphericalTensorN<Cmpt, length>& st
+)
 {
     const Cmpt& s = st.v_[0];
-    VectorSpaceOps<DiagTensorN<Cmpt, length>::nComponents,0>::eqOpS(*this, s, eqOp<Cmpt>());
+    VectorSpaceOps<DiagTensorN<Cmpt, length>::nComponents,0>::eqOpS
+    (
+        *this,
+        s,
+        eqOp<Cmpt>()
+    );
 }
 
 // * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
@@ -105,10 +119,21 @@ inline void DiagTensorN<Cmpt, length>::operator=(const SphericalTensorN<Cmpt, le
 //- Addition of DiagTensorN and DiagTensorN
 template <class Cmpt, int length>
 inline DiagTensorN<Cmpt, length>
-operator+(const DiagTensorN<Cmpt,length>& dt1, const DiagTensorN<Cmpt,length>& dt2)
+operator+
+(
+    const DiagTensorN<Cmpt, length>& dt1,
+    const DiagTensorN<Cmpt, length>& dt2
+)
 {
     DiagTensorN<Cmpt, length> res;
-    VectorSpaceOps<DiagTensorN<Cmpt,length>::nComponents,0>::op(res, dt1, dt2, plusOp<Cmpt>());
+    VectorSpaceOps<DiagTensorN<Cmpt, length>::nComponents,0>::op
+    (
+        res,
+        dt1,
+        dt2,
+        plusOp<Cmpt>()
+    );
+
     return res;
 }
 
@@ -116,11 +141,22 @@ operator+(const DiagTensorN<Cmpt,length>& dt1, const DiagTensorN<Cmpt,length>& d
 //- Addition of DiagTensorN and SphericalTensorN
 template <class Cmpt, int length>
 inline DiagTensorN<Cmpt, length>
-operator+(const DiagTensorN<Cmpt,length>& dt1, const SphericalTensorN<Cmpt,length>& st2)
+operator+
+(
+    const DiagTensorN<Cmpt, length>& dt1,
+    const SphericalTensorN<Cmpt, length>& st2
+)
 {
     const Cmpt& s = st2.v_[0];
     DiagTensorN<Cmpt, length> res;
-    VectorSpaceOps<DiagTensorN<Cmpt,length>::nComponents,0>::opVS(res, dt1, s, plusOp<Cmpt>());
+    VectorSpaceOps<DiagTensorN<Cmpt, length>::nComponents,0>::opVS
+    (
+        res,
+        dt1,
+        s,
+        plusOp<Cmpt>()
+    );
+
     return res;
 }
 
@@ -128,11 +164,22 @@ operator+(const DiagTensorN<Cmpt,length>& dt1, const SphericalTensorN<Cmpt,lengt
 //- Addition of SphericalTensorN and DiagTensorN
 template <class Cmpt, int length>
 inline DiagTensorN<Cmpt, length>
-operator+(const SphericalTensorN<Cmpt,length>& st1, const DiagTensorN<Cmpt,length>& dt2)
+operator+
+(
+    const SphericalTensorN<Cmpt, length>& st1,
+    const DiagTensorN<Cmpt, length>& dt2
+)
 {
     const Cmpt& s = st1.v_[0];
     DiagTensorN<Cmpt, length> res;
-    VectorSpaceOps<DiagTensorN<Cmpt,length>::nComponents,0>::opSV(res, s, dt2, plusOp<Cmpt>());
+    VectorSpaceOps<DiagTensorN<Cmpt, length>::nComponents,0>::opSV
+    (
+        res,
+        s,
+        dt2,
+        plusOp<Cmpt>()
+    );
+
     return res;
 }
 
@@ -140,10 +187,21 @@ operator+(const SphericalTensorN<Cmpt,length>& st1, const DiagTensorN<Cmpt,lengt
 //- Subtraction of DiagTensorN and DiagTensorN
 template <class Cmpt, int length>
 inline DiagTensorN<Cmpt, length>
-operator-(const DiagTensorN<Cmpt,length>& dt1, const DiagTensorN<Cmpt,length>& dt2)
+operator-
+(
+    const DiagTensorN<Cmpt, length>& dt1,
+    const DiagTensorN<Cmpt, length>& dt2
+)
 {
     DiagTensorN<Cmpt, length> res;
-    VectorSpaceOps<DiagTensorN<Cmpt,length>::nComponents,0>::op(res, dt1, dt2, minusOp<Cmpt>());
+    VectorSpaceOps<DiagTensorN<Cmpt, length>::nComponents,0>::op
+    (
+        res,
+        dt1,
+        dt2,
+        minusOp<Cmpt>()
+    );
+
     return res;
 }
 
@@ -151,11 +209,22 @@ operator-(const DiagTensorN<Cmpt,length>& dt1, const DiagTensorN<Cmpt,length>& d
 //- Subtraction of DiagTensorN and SphericalTensorN
 template <class Cmpt, int length>
 inline DiagTensorN<Cmpt, length>
-operator-(const DiagTensorN<Cmpt,length>& dt1, const SphericalTensorN<Cmpt,length>& st2)
+operator-
+(
+    const DiagTensorN<Cmpt, length>& dt1,
+    const SphericalTensorN<Cmpt, length>& st2
+)
 {
     const Cmpt& s = st2.v_[0];
     DiagTensorN<Cmpt, length> res;
-    VectorSpaceOps<DiagTensorN<Cmpt,length>::nComponents,0>::opVS(res, dt1, s, minusOp<Cmpt>());
+    VectorSpaceOps<DiagTensorN<Cmpt, length>::nComponents,0>::opVS
+    (
+        res,
+        dt1,
+        s,
+        minusOp<Cmpt>()
+    );
+
     return res;
 }
 
@@ -163,11 +232,22 @@ operator-(const DiagTensorN<Cmpt,length>& dt1, const SphericalTensorN<Cmpt,lengt
 //- Subtraction of SphericalTensorN and DiagTensorN
 template <class Cmpt, int length>
 inline DiagTensorN<Cmpt, length>
-operator-(const SphericalTensorN<Cmpt,length>& st1, const DiagTensorN<Cmpt,length>& dt2)
+operator-
+(
+    const SphericalTensorN<Cmpt, length>& st1,
+    const DiagTensorN<Cmpt, length>& dt2
+)
 {
     const Cmpt& s = st1.v_[0];
     DiagTensorN<Cmpt, length> res;
-    VectorSpaceOps<DiagTensorN<Cmpt,length>::nComponents,0>::opSV(res, s, dt2, minusOp<Cmpt>());
+    VectorSpaceOps<DiagTensorN<Cmpt, length>::nComponents,0>::opSV
+    (
+        res,
+        s,
+        dt2,
+        minusOp<Cmpt>()
+    );
+
     return res;
 }
 
@@ -176,10 +256,21 @@ operator-(const SphericalTensorN<Cmpt,length>& st1, const DiagTensorN<Cmpt,lengt
 template <class Cmpt, int length>
 inline typename
 innerProduct<DiagTensorN<Cmpt, length>, DiagTensorN<Cmpt, length> >::type
-operator&(const DiagTensorN<Cmpt, length>& dt1, const DiagTensorN<Cmpt, length>& dt2)
+operator&
+(
+    const DiagTensorN<Cmpt, length>& dt1,
+    const DiagTensorN<Cmpt, length>& dt2
+)
 {
     DiagTensorN<Cmpt, length> res;
-    VectorSpaceOps<VectorN<Cmpt, length>::nComponents,0>::op(res, dt1, dt2, multiplyOp<Cmpt>());
+    VectorSpaceOps<VectorN<Cmpt, length>::nComponents,0>::op
+    (
+        res,
+        dt1,
+        dt2,
+        multiplyOp<Cmpt>()
+    );
+
     return res;
 }
 
@@ -188,11 +279,22 @@ operator&(const DiagTensorN<Cmpt, length>& dt1, const DiagTensorN<Cmpt, length>&
 template <class Cmpt, int length>
 inline typename
 innerProduct<SphericalTensorN<Cmpt, length>, DiagTensorN<Cmpt, length> >::type
-operator&(const SphericalTensorN<Cmpt, length>& st1, const DiagTensorN<Cmpt, length>& dt2)
+operator&
+(
+    const SphericalTensorN<Cmpt, length>& st1,
+    const DiagTensorN<Cmpt, length>& dt2
+)
 {
     const Cmpt& s = st1.v_[0];
     DiagTensorN<Cmpt, length> res;
-    VectorSpaceOps<VectorN<Cmpt, length>::nComponents,0>::opSV(res, s, dt2, multiplyOp<Cmpt>());
+    VectorSpaceOps<VectorN<Cmpt, length>::nComponents,0>::opSV
+    (
+        res,
+        s,
+        dt2,
+        multiplyOp<Cmpt>()
+    );
+
     return res;
 }
 
@@ -201,11 +303,22 @@ operator&(const SphericalTensorN<Cmpt, length>& st1, const DiagTensorN<Cmpt, len
 template <class Cmpt, int length>
 inline typename
 innerProduct<DiagTensorN<Cmpt, length>, SphericalTensorN<Cmpt, length> >::type
-operator&(const DiagTensorN<Cmpt, length>& dt1, const SphericalTensorN<Cmpt, length>& st2)
+operator&
+(
+    const DiagTensorN<Cmpt, length>& dt1,
+    const SphericalTensorN<Cmpt, length>& st2
+)
 {
     const Cmpt& s = st2.v_[0];
     DiagTensorN<Cmpt, length> res;
-    VectorSpaceOps<VectorN<Cmpt, length>::nComponents,0>::opVS(res, dt1, s, multiplyOp<Cmpt>());
+    VectorSpaceOps<VectorN<Cmpt, length>::nComponents,0>::opVS
+    (
+        res,
+        dt1,
+        s,
+        multiplyOp<Cmpt>()
+    );
+
     return res;
 }
 
@@ -214,10 +327,21 @@ operator&(const DiagTensorN<Cmpt, length>& dt1, const SphericalTensorN<Cmpt, len
 template <class Cmpt, int length>
 inline typename
 innerProduct<DiagTensorN<Cmpt, length>, VectorN<Cmpt, length> >::type
-operator&(const DiagTensorN<Cmpt, length>& dt, const VectorN<Cmpt, length>& v)
+operator&
+(
+    const DiagTensorN<Cmpt, length>& dt,
+    const VectorN<Cmpt, length>& v
+)
 {
     VectorN<Cmpt, length> res;
-    VectorSpaceOps<VectorN<Cmpt, length>::nComponents,0>::opVV(res, dt, v, multiplyOp<Cmpt>());
+    VectorSpaceOps<VectorN<Cmpt, length>::nComponents,0>::opVV
+    (
+        res,
+        dt,
+        v,
+        multiplyOp<Cmpt>()
+    );
+
     return res;
 }
 
@@ -227,10 +351,21 @@ operator&(const DiagTensorN<Cmpt, length>& dt, const VectorN<Cmpt, length>& v)
 template <class Cmpt, int length>
 inline typename
 innerProduct<VectorN<Cmpt, length>, DiagTensorN<Cmpt, length> >::type
-operator&(const VectorN<Cmpt, length>& v, const DiagTensorN<Cmpt, length>& dt)
+operator&
+(
+    const VectorN<Cmpt, length>& v,
+    const DiagTensorN<Cmpt, length>& dt
+)
 {
     VectorN<Cmpt, length> res;
-    VectorSpaceOps<VectorN<Cmpt, length>::nComponents,0>::opVV(res, v, dt, multiplyOp<Cmpt>());
+    VectorSpaceOps<VectorN<Cmpt, length>::nComponents,0>::opVV
+    (
+        res,
+        v,
+        dt,
+        multiplyOp<Cmpt>()
+    );
+
     return res;
 }
 
@@ -241,7 +376,14 @@ inline DiagTensorN<Cmpt, length>
 operator/(const scalar s, const DiagTensorN<Cmpt, length>& dt)
 {
     DiagTensorN<Cmpt, length> res;
-    VectorSpaceOps<DiagTensorN<Cmpt, length>::nComponents,0>::opSV(res, s, dt, divideOp3<Cmpt, scalar, Cmpt>());
+    VectorSpaceOps<DiagTensorN<Cmpt, length>::nComponents,0>::opSV
+    (
+        res,
+        s,
+        dt,
+        divideOp3<Cmpt, scalar, Cmpt>()
+    );
+
     return res;
 }
 
@@ -252,7 +394,13 @@ inline VectorN<Cmpt,length>
 operator/(const VectorN<Cmpt, length>& v, const DiagTensorN<Cmpt, length>& dt)
 {
     VectorN<Cmpt, length> res(v);
-    VectorSpaceOps<DiagTensorN<Cmpt, length>::nComponents,0>::eqOp(res, dt, divideEqOp<Cmpt>());
+    VectorSpaceOps<DiagTensorN<Cmpt, length>::nComponents,0>::eqOp
+    (
+        res,
+        dt,
+        divideEqOp<Cmpt>()
+    );
+
     return res;
 }
 
@@ -260,10 +408,21 @@ operator/(const VectorN<Cmpt,length>& v, const DiagTensorN<Cmpt,length>& dt)
 //- Inner Product of a DiagTensorN and an inverse DiagTensorN
 template <class Cmpt, int length>
 inline DiagTensorN<Cmpt, length>
-operator/(const DiagTensorN<Cmpt,length>& dt1, const DiagTensorN<Cmpt,length>& dt2)
+operator/
+(
+    const DiagTensorN<Cmpt, length>& dt1,
+    const DiagTensorN<Cmpt, length>& dt2
+)
 {
     DiagTensorN<Cmpt, length> res;
-    VectorSpaceOps<DiagTensorN<Cmpt, length>::nComponents,0>::op(res, dt1, dt2, divideOp<Cmpt>());
+    VectorSpaceOps<DiagTensorN<Cmpt, length>::nComponents,0>::op
+    (
+        res,
+        dt1,
+        dt2,
+        divideOp<Cmpt>()
+    );
+
     return res;
 }
 
@@ -271,11 +430,22 @@ operator/(const DiagTensorN<Cmpt,length>& dt1, const DiagTensorN<Cmpt,length>& d
 //- Inner Product of a SphericalTensorN and an inverse DiagTensorN
 template <class Cmpt, int length>
 inline DiagTensorN<Cmpt, length>
-operator/(const SphericalTensorN<Cmpt,length>& st1, const DiagTensorN<Cmpt,length>& dt2)
+operator/
+(
+    const SphericalTensorN<Cmpt, length>& st1,
+    const DiagTensorN<Cmpt, length>& dt2
+)
 {
     const Cmpt& s = st1.v_[0];
     DiagTensorN<Cmpt, length> res;
-    VectorSpaceOps<DiagTensorN<Cmpt, length>::nComponents,0>::opSV(res, s, dt2, divideOp<Cmpt>());
+    VectorSpaceOps<DiagTensorN<Cmpt, length>::nComponents,0>::opSV
+    (
+        res,
+        s,
+        dt2,
+        divideOp<Cmpt>()
+    );
+
     return res;
 }
 
@@ -283,11 +453,22 @@ operator/(const SphericalTensorN<Cmpt,length>& st1, const DiagTensorN<Cmpt,lengt
 //- Inner Product of a DiagTensorN and an inverse SphericalTensorN
 template <class Cmpt, int length>
 inline DiagTensorN<Cmpt, length>
-operator/(const DiagTensorN<Cmpt,length>& dt1, const SphericalTensorN<Cmpt,length>& st2)
+operator/
+(
+    const DiagTensorN<Cmpt, length>& dt1,
+    const SphericalTensorN<Cmpt, length>& st2
+)
 {
     const Cmpt& s = st2.v_[0];
     DiagTensorN<Cmpt, length> res;
-    VectorSpaceOps<DiagTensorN<Cmpt, length>::nComponents,0>::opVS(res, dt1, s, divideOp<Cmpt>());
+    VectorSpaceOps<DiagTensorN<Cmpt, length>::nComponents,0>::opVS
+    (
+        res,
+        dt1,
+        s,
+        divideOp<Cmpt>()
+    );
+
     return res;
 }
 
@@ -297,7 +478,14 @@ template <class Cmpt, int length>
 inline DiagTensorN<Cmpt, length> inv(const DiagTensorN<Cmpt, length>& dt)
 {
     DiagTensorN<Cmpt, length> res;
-    VectorSpaceOps<DiagTensorN<Cmpt, length>::nComponents,0>::opSV(res, 1.0, dt, divideOp<Cmpt>());
+    VectorSpaceOps<DiagTensorN<Cmpt, length>::nComponents,0>::opSV
+    (
+        res,
+        1.0,
+        dt,
+        divideOp<Cmpt>()
+    );
+
     return res;
 }
 

src/foam/primitives/VectorN/SphericalTensorNI.H

@@ -79,14 +79,16 @@ inline SphericalTensorN<Cmpt, length>::SphericalTensorN(Istream& is)
 
 //- Return diagonal tensor diagonal
 template <class Cmpt, int length>
-inline SphericalTensorN<Cmpt, length> SphericalTensorN<Cmpt, length>::diag() const
+inline SphericalTensorN<Cmpt, length>
+SphericalTensorN<Cmpt, length>::diag() const
 {
     return *this;
 }
 
 //- Return spherical tensor transpose
 template <class Cmpt, int length>
-inline SphericalTensorN<Cmpt, length> SphericalTensorN<Cmpt, length>::T() const
+inline SphericalTensorN<Cmpt, length>
+SphericalTensorN<Cmpt, length>::T() const
 {
     return *this;
 }
@@ -112,7 +114,11 @@ inline SphericalTensorN<Cmpt, length> transform
 //- Addition of SphericalTensorN and SphericalTensorN
 template <class Cmpt, int length>
 inline SphericalTensorN<Cmpt, length>
-operator+(const SphericalTensorN<Cmpt,length>& st1, const SphericalTensorN<Cmpt,length>& st2)
+operator+
+(
+    const SphericalTensorN<Cmpt, length>& st1,
+    const SphericalTensorN<Cmpt, length>& st2
+)
 {
     return SphericalTensorN<Cmpt, length>(st1.v_[0] + st2.v_[0]);
 }
@@ -121,7 +127,11 @@ operator+(const SphericalTensorN<Cmpt,length>& st1, const SphericalTensorN<Cmpt,
 //- Subtraction of SphericalTensorN and SphericalTensorN
 template <class Cmpt, int length>
 inline SphericalTensorN<Cmpt, length>
-operator-(const SphericalTensorN<Cmpt,length>& st1, const SphericalTensorN<Cmpt,length>& st2)
+operator-
+(
+    const SphericalTensorN<Cmpt, length>& st1,
+    const SphericalTensorN<Cmpt, length>& st2
+)
 {
     return SphericalTensorN<Cmpt, length>(st1.v_[0] - st2.v_[0]);
 }
@@ -130,8 +140,16 @@ operator-(const SphericalTensorN<Cmpt,length>& st1, const SphericalTensorN<Cmpt,
 //- Inner-product between spherical tensor and spherical tensor
 template <class Cmpt, int length>
 inline typename
-innerProduct<SphericalTensorN<Cmpt, length>, SphericalTensorN<Cmpt, length> >::type
-operator&(const SphericalTensorN<Cmpt, length>& st1, const SphericalTensorN<Cmpt, length>& st2)
+innerProduct
+<
+    SphericalTensorN<Cmpt, length>,
+    SphericalTensorN<Cmpt, length>
+>::type
+operator&
+(
+    const SphericalTensorN<Cmpt, length>& st1,
+    const SphericalTensorN<Cmpt, length>& st2
+)
 {
     return SphericalTensorN<Cmpt, length>(st1.v_[0]*st2.v_[0]);
 }
@@ -141,11 +159,22 @@ operator&(const SphericalTensorN<Cmpt, length>& st1, const SphericalTensorN<Cmpt
 template <class Cmpt, int length>
 inline typename
 innerProduct<SphericalTensorN<Cmpt, length>, VectorN<Cmpt, length> >::type
-operator&(const SphericalTensorN<Cmpt, length>& st, const VectorN<Cmpt, length>& v)
+operator&
+(
+    const SphericalTensorN<Cmpt, length>& st,
+    const VectorN<Cmpt, length>& v
+)
 {
     const Cmpt& s = st.v_[0];
     VectorN<Cmpt, length> res;
-    VectorSpaceOps<VectorN<Cmpt, length>::nComponents,0>::opSV(res, s, v, multiplyOp<Cmpt>());
+    VectorSpaceOps<VectorN<Cmpt, length>::nComponents,0>::opSV
+    (
+        res,
+        s,
+        v,
+        multiplyOp<Cmpt>()
+    );
+
     return res;
 }
 
@@ -154,11 +183,22 @@ operator&(const SphericalTensorN<Cmpt, length>& st, const VectorN<Cmpt, length>&
 template <class Cmpt, int length>
 inline typename
 innerProduct<VectorN<Cmpt, length>, SphericalTensorN<Cmpt, length> >::type
-operator&(const VectorN<Cmpt, length>& v, const SphericalTensorN<Cmpt, length>& st)
+operator&
+(
+    const VectorN<Cmpt, length>& v,
+    const SphericalTensorN<Cmpt, length>& st
+)
 {
     const Cmpt& s = st.v_[0];
     VectorN<Cmpt, length> res;
-    VectorSpaceOps<VectorN<Cmpt, length>::nComponents,0>::opVS(res, v, s, multiplyOp<Cmpt>());
+    VectorSpaceOps<VectorN<Cmpt, length>::nComponents,0>::opVS
+    (
+        res,
+        v,
+        s,
+        multiplyOp<Cmpt>()
+    );
+
     return res;
 }
 
@@ -184,7 +224,11 @@ operator/(const scalar s, const SphericalTensorN<Cmpt, length>& st)
 //- Inner Product of a VectorN by an inverse SphericalTensorN
 template <class Cmpt, int length>
 inline VectorN<Cmpt, length>
-operator/(const VectorN<Cmpt,length>& v, const SphericalTensorN<Cmpt,length>& st)
+operator/
+(
+    const VectorN<Cmpt, length>& v,
+    const SphericalTensorN<Cmpt, length>& st
+)
 {
     return v/st.v_[0];
 }
@@ -193,7 +237,11 @@ operator/(const VectorN<Cmpt,length>& v, const SphericalTensorN<Cmpt,length>& st
 //- Inner Product of a SphericalTensorN and an inverse SphericalTensorN
 template <class Cmpt, int length>
 inline SphericalTensorN<Cmpt, length>
-operator/(const SphericalTensorN<Cmpt,length>& st1, const SphericalTensorN<Cmpt,length>& st2)
+operator/
+(
+    const SphericalTensorN<Cmpt, length>& st1,
+    const SphericalTensorN<Cmpt, length>& st2
+)
 {
     return SphericalTensorN<Cmpt, length>(st1.v_[0]/st2.v_[0]);
 }
@@ -201,7 +249,10 @@ operator/(const SphericalTensorN<Cmpt,length>& st1, const SphericalTensorN<Cmpt,
 
 //- Return the inverse of a spherical tensor
 template <class Cmpt, int length>
-inline SphericalTensorN<Cmpt, length> inv(const SphericalTensorN<Cmpt, length>& st)
+inline SphericalTensorN<Cmpt, length> inv
+(
+    const SphericalTensorN<Cmpt, length>& st
+)
 {
     return SphericalTensorN<Cmpt, length>(pTraits<Cmpt>::one/st.v_[0]);
 }
@@ -209,7 +260,10 @@ inline SphericalTensorN<Cmpt, length> inv(const SphericalTensorN<Cmpt, length>&
 
 //- Return tensor diagonal
 template <class Cmpt, int length>
-inline SphericalTensorN<Cmpt, length> diag(const SphericalTensorN<Cmpt, length>& st)
+inline SphericalTensorN<Cmpt, length> diag
+(
+    const SphericalTensorN<Cmpt, length>& st
+)
 {
     return st;
 }
@@ -258,7 +312,11 @@ public:
 };
 
 template<class Cmpt, int length>
-class innerProduct<SphericalTensorN<Cmpt, length>, SphericalTensorN<Cmpt, length> >
+class innerProduct
+<
+    SphericalTensorN<Cmpt, length>,
+    SphericalTensorN<Cmpt, length>
+>
 {
 public: