Third version of BlockILUCp preconditioner - still not working properly

src/foam/Make/files

@@ -663,7 +663,11 @@ matrices/blockLduMatrix/BlockLduMatrix/symmTensorBlockLduMatrix.C
 matrices/blockLduMatrix/BlockLduMatrix/tensorBlockLduMatrix.C
 matrices/blockLduMatrix/BlockLduMatrix/blockVectorNMatrices.C
 matrices/blockLduMatrix/BlockLduMatrix/BlockConstraint/scalarBlockConstraint.C
+
 matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrices.C
+matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/sphericalTensorExtendedBlockLduMatrix.C
+matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/symmTensorExtendedBlockLduMatrix.C
+matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/tensorExtendedBlockLduMatrix.C
 matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockVectorNMatrices.C
 
 BlockLduInterfaceField = matrices/blockLduMatrix/BlockLduMatrix/BlockLduInterfaceField

src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrices.H

@@ -35,7 +35,9 @@ Author
 #ifndef extendedBlockLduMatrices_H
 #define extendedBlockLduMatrices_H
 
-#include "extendedBlockLduMatrix.H"
+#include "sphericalTensorExtendedBlockLduMatrix.H"
+#include "symmTensorExtendedBlockLduMatrix.H"
+#include "tensorExtendedBlockLduMatrix.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.C

@@ -25,36 +25,14 @@ License
 
 #include "extendedBlockLduMatrix.H"
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<class Type>
-Foam::extendedBlockLduMatrix<Type>::extendedBlockLduMatrix
+void Foam::extendedBlockLduMatrix<Type>::mapOffDiagCoeffs
 (
-    const BlockLduMatrix<Type>& blockLdum,
+    const BlockLduMatrix<Type>& blockLdum
-    const label extensionLevel,
-    const polyMesh& polyMesh
 )
-:
-    basicBlockLduMatrix_(blockLdum),
-    extLduAddr_
-    (
-        extendedLduAddressing::New
-        (
-            polyMesh,
-            blockLdum.lduAddr(),
-            extensionLevel
-        )
-    ),
-    extendedLowerPtr_(NULL),
-    extendedUpperPtr_(NULL)
 {
-    if (debug)
-    {
-        Info<< "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&) :"
-               "Constructing extendedBlockLduMatrix."
-            << endl;
-    }
-
     if (blockLdum.diagonal())
     {
         WarningIn("extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)")
@@ -299,6 +277,42 @@ Foam::extendedBlockLduMatrix<Type>::extendedBlockLduMatrix
 }
 
 
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::extendedBlockLduMatrix<Type>::extendedBlockLduMatrix
+(
+    const BlockLduMatrix<Type>& blockLdum,
+    const label extensionLevel,
+    const polyMesh& polyMesh
+)
+:
+    basicBlockLduMatrix_(blockLdum),
+    extLduAddr_
+    (
+        extendedLduAddressing::New
+        (
+            polyMesh,
+            blockLdum.lduAddr(),
+            extensionLevel
+        )
+    ),
+    extendedLowerPtr_(NULL),
+    extendedUpperPtr_(NULL)
+{
+    if (debug)
+    {
+        Info<< "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&) :"
+               "Constructing extendedBlockLduMatrix."
+            << endl;
+    }
+
+    // Map off diag coeffs from original block matrix to this extended block
+    // matrix
+    mapOffDiagCoeffs(blockLdum);
+}
+
+
 // * * * * * * * * * * * * * * * * Destructor* * * * * * * * * * * * * * * * //
 
 template<class Type>

src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.H

@@ -86,6 +86,10 @@ private:
         //- Disallow default bitwise assignement
         void operator=(const extendedBlockLduMatrix<Type>&);
 
+        //- Map upper and lower coeffs from ordinary block matrix to extended
+        //  block matrix
+        void mapOffDiagCoeffs(const BlockLduMatrix<Type>&);
+
 
 public:
 

src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/sphericalTensorExtendedBlockLduMatrix.C

@@ -0,0 +1,249 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     |
+    \\  /    A nd           | For copyright notice see file Copyright
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 sphericalTensorExtendedBlockLduMatrix_H
+#define sphericalTensorExtendedBlockLduMatrix_H
+
+#include "extendedBlockLduMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<>
+void Foam::extendedBlockLduMatrix<Foam::sphericalTensor>::mapOffDiagCoeffs
+(
+    const BlockLduMatrix<sphericalTensor>& blockLdum
+)
+{
+    if (blockLdum.diagonal())
+    {
+        WarningIn("extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)")
+            << "Attempted to create extended lower/upper coeffs for block "
+            << "matrix that is diagonal."
+            << nl << endl;
+    }
+    else if (blockLdum.symmetric())
+    {
+        // Get reference to faceMap in extended addressing
+        const unallocLabelList& faceMap = extLduAddr_.faceMap();
+
+        // Avoid assuming it's upper if the matrix is symmetric
+        if (blockLdum.thereIsUpper())
+        {
+            // Allocate extended upper only
+            extendedUpperPtr_ = new TypeCoeffField
+            (
+                extLduAddr_.extendedUpperAddr().size()
+            );
+            TypeCoeffField& extUpper = *extendedUpperPtr_;
+
+            // Get upper coeffs from underlying lduMatrix
+            const TypeCoeffField& upper = blockLdum.upper();
+
+            if (upper.activeType() == blockCoeffBase::SCALAR)
+            {
+                // Helper type definition
+                typedef typename CoeffField<sphericalTensor>::scalarTypeField
+                    activeType;
+
+                // Get references to fields
+                const activeType& activeUpper = upper.asScalar();
+                activeType& activeExtUpper = extUpper.asScalar();
+
+                // Copy non-zero coeffs from basic lduMatrix into corresponding
+                // positions
+                forAll (upper, faceI)
+                {
+                    activeExtUpper[faceMap[faceI]] = activeUpper[faceI];
+                }
+            }
+            else if (upper.activeType() == blockCoeffBase::LINEAR)
+            {
+                // Helper type definition
+                typedef typename CoeffField<sphericalTensor>::linearTypeField
+                    activeType;
+
+                // Get references to fields
+                const activeType& activeUpper = upper.asLinear();
+                activeType& activeExtUpper = extUpper.asLinear();
+
+                // Copy non-zero coeffs from basic lduMatrix into corresponding
+                // positions
+                forAll (upper, faceI)
+                {
+                    activeExtUpper[faceMap[faceI]] = activeUpper[faceI];
+                }
+            }
+            else
+            {
+                FatalErrorIn
+                (
+                    "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
+                )   << "Problem between ordinary block matrix and extended"
+                    << " block matrix upper coeffs type morphing."
+                    << abort(FatalError);
+            }
+        }
+        else
+        {
+            // Allocate extended lower only
+            extendedLowerPtr_ = new TypeCoeffField
+            (
+                extLduAddr_.extendedLowerAddr().size()
+            );
+            TypeCoeffField& extLower = *extendedLowerPtr_;
+
+            // Get lower coeffs from underlying lduMatrix
+            const TypeCoeffField& lower = blockLdum.lower();
+
+            if (lower.activeType() == blockCoeffBase::SCALAR)
+            {
+                // Helper type definition
+                typedef typename CoeffField<sphericalTensor>::scalarTypeField
+                    activeType;
+
+                // Get references to fields
+                const activeType& activeLower = lower.asScalar();
+                activeType& activeExtLower = extLower.asScalar();
+
+                // Copy non-zero coeffs from basic lduMatrix into corresponding
+                // positions
+                forAll (lower, faceI)
+                {
+                    activeExtLower[faceMap[faceI]] = activeLower[faceI];
+                }
+            }
+            else if (lower.activeType() == blockCoeffBase::LINEAR)
+            {
+                // Helper type definition
+                typedef typename CoeffField<sphericalTensor>::linearTypeField
+                    activeType;
+
+                // Get references to fields
+                const activeType& activeLower = lower.asLinear();
+                activeType& activeExtLower = extLower.asLinear();
+
+                // Copy non-zero coeffs from basic lduMatrix into corresponding
+                // positions
+                forAll (lower, faceI)
+                {
+                    activeExtLower[faceMap[faceI]] = activeLower[faceI];
+                }
+            }
+            else
+            {
+                FatalErrorIn
+                (
+                    "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
+                )   << "Problem between ordinary block matrix and extended"
+                    << " block matrix lower coeffs type morphing."
+                    << abort(FatalError);
+            }
+        }
+    }
+    else
+    {
+        // Get reference to faceMap in extended addressing
+        const unallocLabelList& faceMap = extLduAddr_.faceMap();
+
+        // Get number of extended faces
+        const label nExtFaces = extLduAddr_.extendedUpperAddr().size();
+
+        // Allocate extended upper and lower
+        extendedUpperPtr_ = new TypeCoeffField(nExtFaces);
+        TypeCoeffField& extUpper = *extendedUpperPtr_;
+
+        extendedLowerPtr_ = new TypeCoeffField(nExtFaces);
+        TypeCoeffField& extLower = *extendedLowerPtr_;
+
+        // Get upper and lower coeffs from underlying lduMatrix
+        const TypeCoeffField& upper = blockLdum.upper();
+        const TypeCoeffField& lower = blockLdum.lower();
+
+        // Assuming lower and upper have the same active type
+        if (upper.activeType() == blockCoeffBase::SCALAR)
+        {
+            // Helper type definition
+            typedef typename CoeffField<sphericalTensor>::scalarTypeField
+                activeType;
+
+            // Get references to fields
+            const activeType& activeUpper = upper.asScalar();
+            activeType& activeExtUpper = extUpper.asScalar();
+            const activeType& activeLower = lower.asScalar();
+            activeType& activeExtLower = extLower.asScalar();
+
+            // Copy non-zero coeffs from basic lduMatrix into corresponding
+            // positions
+            forAll (upper, faceI)
+            {
+                activeExtUpper[faceMap[faceI]] = activeUpper[faceI];
+                activeExtLower[faceMap[faceI]] = activeLower[faceI];
+            }
+        }
+        else if (upper.activeType() == blockCoeffBase::LINEAR)
+        {
+            // Helper type definition
+            typedef typename CoeffField<sphericalTensor>::linearTypeField
+                activeType;
+
+            // Get references to fields
+            const activeType& activeUpper = upper.asLinear();
+            activeType& activeExtUpper = extUpper.asLinear();
+            const activeType& activeLower = lower.asLinear();
+            activeType& activeExtLower = extLower.asLinear();
+
+            // Copy non-zero coeffs from basic lduMatrix into corresponding
+            // positions
+            forAll (upper, faceI)
+            {
+                activeExtUpper[faceMap[faceI]] = activeUpper[faceI];
+                activeExtLower[faceMap[faceI]] = activeLower[faceI];
+            }
+        }
+        else
+        {
+            FatalErrorIn
+            (
+                "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
+            )   << "Problem between ordinary block matrix and extended"
+                << " block matrix upper/lower coeffs type morphing."
+                << abort(FatalError);
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/sphericalTensorExtendedBlockLduMatrix.H

@@ -0,0 +1,63 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     |
+    \\  /    A nd           | For copyright notice see file Copyright
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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
+    extendedBlockLduMatrix
+
+Description
+    Template specialisation for sphericalTensor extended block matrix
+
+Author
+    Vuko Vukcevic.  All rights reserved.
+
+SourceFiles
+    sphericalTensorExtendedBlockLduMatrix.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef sphericalTensorExtendedBlockLduMatrix_H
+#define sphericalTensorExtendedBlockLduMatrix_H
+
+#include "extendedBlockLduMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<>
+void extendedBlockLduMatrix<sphericalTensor>::mapOffDiagCoeffs
+(
+    const BlockLduMatrix<sphericalTensor>& blockLdum
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/symmTensorExtendedBlockLduMatrix.C

@@ -0,0 +1,247 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     |
+    \\  /    A nd           | For copyright notice see file Copyright
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 symmTensorExtendedBlockLduMatrix_H
+#define symmTensorExtendedBlockLduMatrix_H
+
+#include "extendedBlockLduMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<>
+void Foam::extendedBlockLduMatrix<Foam::symmTensor>::mapOffDiagCoeffs
+(
+    const BlockLduMatrix<symmTensor>& blockLdum
+)
+{
+    if (blockLdum.diagonal())
+    {
+        WarningIn("extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)")
+            << "Attempted to create extended lower/upper coeffs for block "
+            << "matrix that is diagonal."
+            << nl << endl;
+    }
+    else if (blockLdum.symmetric())
+    {
+        // Get reference to faceMap in extended addressing
+        const unallocLabelList& faceMap = extLduAddr_.faceMap();
+
+        // Avoid assuming it's upper if the matrix is symmetric
+        if (blockLdum.thereIsUpper())
+        {
+            // Allocate extended upper only
+            extendedUpperPtr_ = new TypeCoeffField
+            (
+                extLduAddr_.extendedUpperAddr().size()
+            );
+            TypeCoeffField& extUpper = *extendedUpperPtr_;
+
+            // Get upper coeffs from underlying lduMatrix
+            const TypeCoeffField& upper = blockLdum.upper();
+
+            if (upper.activeType() == blockCoeffBase::SCALAR)
+            {
+                // Helper type definition
+                typedef typename CoeffField<symmTensor>::scalarTypeField
+                    activeType;
+
+                // Get references to fields
+                const activeType& activeUpper = upper.asScalar();
+                activeType& activeExtUpper = extUpper.asScalar();
+
+                // Copy non-zero coeffs from basic lduMatrix into corresponding
+                // positions
+                forAll (upper, faceI)
+                {
+                    activeExtUpper[faceMap[faceI]] = activeUpper[faceI];
+                }
+            }
+            else if (upper.activeType() == blockCoeffBase::LINEAR)
+            {
+                // Helper type definition
+                typedef typename CoeffField<symmTensor>::linearTypeField
+                    activeType;
+
+                // Get references to fields
+                const activeType& activeUpper = upper.asLinear();
+                activeType& activeExtUpper = extUpper.asLinear();
+
+                // Copy non-zero coeffs from basic lduMatrix into corresponding
+                // positions
+                forAll (upper, faceI)
+                {
+                    activeExtUpper[faceMap[faceI]] = activeUpper[faceI];
+                }
+            }
+            else
+            {
+                FatalErrorIn
+                (
+                    "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
+                )   << "Problem between ordinary block matrix and extended"
+                    << " block matrix upper coeffs type morphing."
+                    << abort(FatalError);
+            }
+        }
+        else
+        {
+            // Allocate extended lower only
+            extendedLowerPtr_ = new TypeCoeffField
+            (
+                extLduAddr_.extendedLowerAddr().size()
+            );
+            TypeCoeffField& extLower = *extendedLowerPtr_;
+
+            // Get lower coeffs from underlying lduMatrix
+            const TypeCoeffField& lower = blockLdum.lower();
+
+            if (lower.activeType() == blockCoeffBase::SCALAR)
+            {
+                // Helper type definition
+                typedef typename CoeffField<symmTensor>::scalarTypeField
+                    activeType;
+
+                // Get references to fields
+                const activeType& activeLower = lower.asScalar();
+                activeType& activeExtLower = extLower.asScalar();
+
+                // Copy non-zero coeffs from basic lduMatrix into corresponding
+                // positions
+                forAll (lower, faceI)
+                {
+                    activeExtLower[faceMap[faceI]] = activeLower[faceI];
+                }
+            }
+            else if (lower.activeType() == blockCoeffBase::LINEAR)
+            {
+                // Helper type definition
+                typedef typename CoeffField<symmTensor>::linearTypeField
+                    activeType;
+
+                // Get references to fields
+                const activeType& activeLower = lower.asLinear();
+                activeType& activeExtLower = extLower.asLinear();
+
+                // Copy non-zero coeffs from basic lduMatrix into corresponding
+                // positions
+                forAll (lower, faceI)
+                {
+                    activeExtLower[faceMap[faceI]] = activeLower[faceI];
+                }
+            }
+            else
+            {
+                FatalErrorIn
+                (
+                    "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
+                )   << "Problem between ordinary block matrix and extended"
+                    << " block matrix lower coeffs type morphing."
+                    << abort(FatalError);
+            }
+        }
+    }
+    else
+    {
+        // Get reference to faceMap in extended addressing
+        const unallocLabelList& faceMap = extLduAddr_.faceMap();
+
+        // Get number of extended faces
+        const label nExtFaces = extLduAddr_.extendedUpperAddr().size();
+
+        // Allocate extended upper and lower
+        extendedUpperPtr_ = new TypeCoeffField(nExtFaces);
+        TypeCoeffField& extUpper = *extendedUpperPtr_;
+
+        extendedLowerPtr_ = new TypeCoeffField(nExtFaces);
+        TypeCoeffField& extLower = *extendedLowerPtr_;
+
+        // Get upper and lower coeffs from underlying lduMatrix
+        const TypeCoeffField& upper = blockLdum.upper();
+        const TypeCoeffField& lower = blockLdum.lower();
+
+        // Assuming lower and upper have the same active type
+        if (upper.activeType() == blockCoeffBase::SCALAR)
+        {
+            // Helper type definition
+            typedef typename CoeffField<symmTensor>::scalarTypeField activeType;
+
+            // Get references to fields
+            const activeType& activeUpper = upper.asScalar();
+            activeType& activeExtUpper = extUpper.asScalar();
+            const activeType& activeLower = lower.asScalar();
+            activeType& activeExtLower = extLower.asScalar();
+
+            // Copy non-zero coeffs from basic lduMatrix into corresponding
+            // positions
+            forAll (upper, faceI)
+            {
+                activeExtUpper[faceMap[faceI]] = activeUpper[faceI];
+                activeExtLower[faceMap[faceI]] = activeLower[faceI];
+            }
+        }
+        else if (upper.activeType() == blockCoeffBase::LINEAR)
+        {
+            // Helper type definition
+            typedef typename CoeffField<symmTensor>::linearTypeField activeType;
+
+            // Get references to fields
+            const activeType& activeUpper = upper.asLinear();
+            activeType& activeExtUpper = extUpper.asLinear();
+            const activeType& activeLower = lower.asLinear();
+            activeType& activeExtLower = extLower.asLinear();
+
+            // Copy non-zero coeffs from basic lduMatrix into corresponding
+            // positions
+            forAll (upper, faceI)
+            {
+                activeExtUpper[faceMap[faceI]] = activeUpper[faceI];
+                activeExtLower[faceMap[faceI]] = activeLower[faceI];
+            }
+        }
+        else
+        {
+            FatalErrorIn
+            (
+                "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
+            )   << "Problem between ordinary block matrix and extended"
+                << " block matrix upper/lower coeffs type morphing."
+                << abort(FatalError);
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/symmTensorExtendedBlockLduMatrix.H

@@ -0,0 +1,63 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     |
+    \\  /    A nd           | For copyright notice see file Copyright
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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
+    extendedBlockLduMatrix
+
+Description
+    Template specialisation for symmTensor extended block matrix
+
+Author
+    Vuko Vukcevic.  All rights reserved.
+
+SourceFiles
+    symmTensorExtendedBlockLduMatrix.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef symmTensorExtendedBlockLduMatrix_H
+#define symmTensorExtendedBlockLduMatrix_H
+
+#include "extendedBlockLduMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<>
+void extendedBlockLduMatrix<symmTensor>::mapOffDiagCoeffs
+(
+    const BlockLduMatrix<symmTensor>& blockLdum
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/tensorExtendedBlockLduMatrix.C

@@ -0,0 +1,243 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     |
+    \\  /    A nd           | For copyright notice see file Copyright
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 tensorExtendedBlockLduMatrix_H
+#define tensorExtendedBlockLduMatrix_H
+
+#include "extendedBlockLduMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<>
+void Foam::extendedBlockLduMatrix<Foam::tensor>::mapOffDiagCoeffs
+(
+    const BlockLduMatrix<tensor>& blockLdum
+)
+{
+    if (blockLdum.diagonal())
+    {
+        WarningIn("extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)")
+            << "Attempted to create extended lower/upper coeffs for block "
+            << "matrix that is diagonal."
+            << nl << endl;
+    }
+    else if (blockLdum.symmetric())
+    {
+        // Get reference to faceMap in extended addressing
+        const unallocLabelList& faceMap = extLduAddr_.faceMap();
+
+        // Avoid assuming it's upper if the matrix is symmetric
+        if (blockLdum.thereIsUpper())
+        {
+            // Allocate extended upper only
+            extendedUpperPtr_ = new TypeCoeffField
+            (
+                extLduAddr_.extendedUpperAddr().size()
+            );
+            TypeCoeffField& extUpper = *extendedUpperPtr_;
+
+            // Get upper coeffs from underlying lduMatrix
+            const TypeCoeffField& upper = blockLdum.upper();
+
+            if (upper.activeType() == blockCoeffBase::SCALAR)
+            {
+                // Helper type definition
+                typedef typename CoeffField<tensor>::scalarTypeField activeType;
+
+                // Get references to fields
+                const activeType& activeUpper = upper.asScalar();
+                activeType& activeExtUpper = extUpper.asScalar();
+
+                // Copy non-zero coeffs from basic lduMatrix into corresponding
+                // positions
+                forAll (upper, faceI)
+                {
+                    activeExtUpper[faceMap[faceI]] = activeUpper[faceI];
+                }
+            }
+            else if (upper.activeType() == blockCoeffBase::LINEAR)
+            {
+                // Helper type definition
+                typedef typename CoeffField<tensor>::linearTypeField activeType;
+
+                // Get references to fields
+                const activeType& activeUpper = upper.asLinear();
+                activeType& activeExtUpper = extUpper.asLinear();
+
+                // Copy non-zero coeffs from basic lduMatrix into corresponding
+                // positions
+                forAll (upper, faceI)
+                {
+                    activeExtUpper[faceMap[faceI]] = activeUpper[faceI];
+                }
+            }
+            else
+            {
+                FatalErrorIn
+                (
+                    "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
+                )   << "Problem between ordinary block matrix and extended"
+                    << " block matrix upper coeffs type morphing."
+                    << abort(FatalError);
+            }
+        }
+        else
+        {
+            // Allocate extended lower only
+            extendedLowerPtr_ = new TypeCoeffField
+            (
+                extLduAddr_.extendedLowerAddr().size()
+            );
+            TypeCoeffField& extLower = *extendedLowerPtr_;
+
+            // Get lower coeffs from underlying lduMatrix
+            const TypeCoeffField& lower = blockLdum.lower();
+
+            if (lower.activeType() == blockCoeffBase::SCALAR)
+            {
+                // Helper type definition
+                typedef typename CoeffField<tensor>::scalarTypeField activeType;
+
+                // Get references to fields
+                const activeType& activeLower = lower.asScalar();
+                activeType& activeExtLower = extLower.asScalar();
+
+                // Copy non-zero coeffs from basic lduMatrix into corresponding
+                // positions
+                forAll (lower, faceI)
+                {
+                    activeExtLower[faceMap[faceI]] = activeLower[faceI];
+                }
+            }
+            else if (lower.activeType() == blockCoeffBase::LINEAR)
+            {
+                // Helper type definition
+                typedef typename CoeffField<tensor>::linearTypeField activeType;
+
+                // Get references to fields
+                const activeType& activeLower = lower.asLinear();
+                activeType& activeExtLower = extLower.asLinear();
+
+                // Copy non-zero coeffs from basic lduMatrix into corresponding
+                // positions
+                forAll (lower, faceI)
+                {
+                    activeExtLower[faceMap[faceI]] = activeLower[faceI];
+                }
+            }
+            else
+            {
+                FatalErrorIn
+                (
+                    "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
+                )   << "Problem between ordinary block matrix and extended"
+                    << " block matrix lower coeffs type morphing."
+                    << abort(FatalError);
+            }
+        }
+    }
+    else
+    {
+        // Get reference to faceMap in extended addressing
+        const unallocLabelList& faceMap = extLduAddr_.faceMap();
+
+        // Get number of extended faces
+        const label nExtFaces = extLduAddr_.extendedUpperAddr().size();
+
+        // Allocate extended upper and lower
+        extendedUpperPtr_ = new TypeCoeffField(nExtFaces);
+        TypeCoeffField& extUpper = *extendedUpperPtr_;
+
+        extendedLowerPtr_ = new TypeCoeffField(nExtFaces);
+        TypeCoeffField& extLower = *extendedLowerPtr_;
+
+        // Get upper and lower coeffs from underlying lduMatrix
+        const TypeCoeffField& upper = blockLdum.upper();
+        const TypeCoeffField& lower = blockLdum.lower();
+
+        // Assuming lower and upper have the same active type
+        if (upper.activeType() == blockCoeffBase::SCALAR)
+        {
+            // Helper type definition
+            typedef typename CoeffField<tensor>::scalarTypeField activeType;
+
+            // Get references to fields
+            const activeType& activeUpper = upper.asScalar();
+            activeType& activeExtUpper = extUpper.asScalar();
+            const activeType& activeLower = lower.asScalar();
+            activeType& activeExtLower = extLower.asScalar();
+
+            // Copy non-zero coeffs from basic lduMatrix into corresponding
+            // positions
+            forAll (upper, faceI)
+            {
+                activeExtUpper[faceMap[faceI]] = activeUpper[faceI];
+                activeExtLower[faceMap[faceI]] = activeLower[faceI];
+            }
+        }
+        else if (upper.activeType() == blockCoeffBase::LINEAR)
+        {
+            // Helper type definition
+            typedef typename CoeffField<tensor>::linearTypeField activeType;
+
+            // Get references to fields
+            const activeType& activeUpper = upper.asLinear();
+            activeType& activeExtUpper = extUpper.asLinear();
+            const activeType& activeLower = lower.asLinear();
+            activeType& activeExtLower = extLower.asLinear();
+
+            // Copy non-zero coeffs from basic lduMatrix into corresponding
+            // positions
+            forAll (upper, faceI)
+            {
+                activeExtUpper[faceMap[faceI]] = activeUpper[faceI];
+                activeExtLower[faceMap[faceI]] = activeLower[faceI];
+            }
+        }
+        else
+        {
+            FatalErrorIn
+            (
+                "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
+            )   << "Problem between ordinary block matrix and extended"
+                << " block matrix upper/lower coeffs type morphing."
+                << abort(FatalError);
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/tensorExtendedBlockLduMatrix.H

@@ -0,0 +1,63 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     |
+    \\  /    A nd           | For copyright notice see file Copyright
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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
+    extendedBlockLduMatrix
+
+Description
+    Template specialisation for tensor extended block matrix
+
+Author
+    Vuko Vukcevic.  All rights reserved.
+
+SourceFiles
+    tensorExtendedBlockLduMatrix.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef tensorExtendedBlockLduMatrix_H
+#define tensorExtendedBlockLduMatrix_H
+
+#include "extendedBlockLduMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<>
+void extendedBlockLduMatrix<tensor>::mapOffDiagCoeffs
+(
+    const BlockLduMatrix<tensor>& blockLdum
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -41,20 +41,25 @@ Author
 
 template<class Type>
 template<class LDUType>
-void Foam::BlockILUCpPrecon<Type>::calcFactorization
+void Foam::BlockILUCpPrecon<Type>::calcActiveTypeFactorization
 (
     Field<LDUType>& preconD,
     Field<LDUType>& extUpper,
-    Field<LDUType>& extLower,
+    Field<LDUType>& extLower
-    Field<LDUType>& zDiag,
+) const
-    Field<LDUType>& z,
-    Field<LDUType>& w
-)
 {
     if (!this->matrix_.diagonal())
     {
+        // Get number of rows
+        const label nRows = preconD.size();
+
+        // Allocate working fields
+        Field<LDUType> z(nRows, pTraits<LDUType>::zero);
+        Field<LDUType> w(nRows, pTraits<LDUType>::zero);
+        LDUType zDiag(pTraits<LDUType>::zero);
+
         // Create multiplication function object
-        typename BlockCoeff<LDUType>::multiply mult;
+        typename BlockCoeff<Type>::multiply mult;
 
         // Get necessary const access to extended ldu addressing
         const extendedLduAddressing& addr = extBlockMatrix_.extendedLduAddr();
@@ -81,10 +86,6 @@ void Foam::BlockILUCpPrecon<Type>::calcFactorization
         // Get access to working fields
         LDUType* __restrict__ zPtr = z.begin();
         LDUType* __restrict__ wPtr = w.begin();
-        LDUType& zDiagI = zDiag[0];
-
-        // Get number of rows
-        const label nRows = preconD.size();
 
         // Define start and end face ("virtual" face when extended addressing is
         // used) of this row/column.
@@ -98,8 +99,9 @@ void Foam::BlockILUCpPrecon<Type>::calcFactorization
             // Start and end of k-th row (upper) and k-th column (lower)
             fStart = ownStartPtr[rowI];
             fEnd = ownStartPtr[rowI + 1];
+
             // Initialize temporary working diagonal
-            zDiagI = diagPtr[rowI];
+            zDiag = diagPtr[rowI];
 
             // Initialize temporary working row and column fields
             for (register label faceI = fStart; faceI < fEnd; ++faceI)
@@ -132,7 +134,7 @@ void Foam::BlockILUCpPrecon<Type>::calcFactorization
                 // Update diagonal
                 // WARNING: Not sure about order of multiplication.
                 // Check it. VV, 3/Jul/2015
-                zDiagI -= mult
+                zDiag -= mult.activeTypeMultiply
                 (
                     lowerPtr[losortCoeff],
                     upperPtr[losortCoeff]
@@ -153,7 +155,7 @@ void Foam::BlockILUCpPrecon<Type>::calcFactorization
                 {
                     // WARNING: Not sure about order of multiplication.
                     // Check it. VV, 3/Jul/2015
-                    zPtr[uPtr[faceI]] -= mult
+                    zPtr[uPtr[faceI]] -= mult.activeTypeMultiply
                     (
                         lowerPtr[losortCoeff],
                         upperPtr[faceI]
@@ -161,17 +163,22 @@ void Foam::BlockILUCpPrecon<Type>::calcFactorization
 
                     // WARNING: Not sure about order of multiplication.
                     // Check it. VV, 3/Jul/2015
-                    wPtr[uPtr[faceI]] -= mult
+                    wPtr[uPtr[faceI]] -= mult.activeTypeMultiply
                     (
-                        upperPtr[losortCoeff],
+                        lowerPtr[faceI],
-                        lowerPtr[faceI]
+                        mult.activeTypeMultiply
+                        (
+//                            diagPtr[lPtr[losortCoeff]], // Note: diag already inverted
+                            mult.inverse(zDiag), // Note: diag already inverted
+                            upperPtr[losortCoeff]
+                        )
                     );
                 }
             }
 
             // Update diagonal entry, inverting it for future use
             LDUType& diagRowI = diagPtr[rowI];
-            diagRowI = mult.inverse(zDiagI);
+            diagRowI = mult.inverse(zDiag);
 
             // Index for updating L and U
             register label zwIndex;
@@ -187,15 +194,16 @@ void Foam::BlockILUCpPrecon<Type>::calcFactorization
 
                 // WARNING: Not sure about order of multiplication.
                 // Check it. VV, 3/Jul/2015
-                lowerPtr[faceI] = mult
+//                lowerPtr[faceI] = mult.activeTypeMultiply
-                (
+//                (
-                    diagRowI,
+//                    diagRowI,
-                    wPtr[zwIndex]
+//                    wPtr[zwIndex]
-                );
+//                );
+                lowerPtr[faceI] = wPtr[zwIndex];
             }
 
             // Reset temporary working fields
-            zDiagI = pTraits<LDUType>::zero;
+            zDiag = pTraits<LDUType>::zero;
 
             // Only reset parts of the working fields that have been updated in
             // this step (for this row and column)
@@ -234,7 +242,7 @@ void Foam::BlockILUCpPrecon<Type>::calcFactorization
         (
             "template <class Type>\n"
             "template <class LDUType>\n"
-            "void BlockILUCpPrecon<Type>::calcFactorization\n"
+            "void BlockILUCpPrecon<Type>::calcActiveTypeFactorization\n"
             "(\n"
             "    Field<LDUType>& preconD,\n"
             "    Field<LDUType>& extUpper,\n"
@@ -252,6 +260,109 @@ void Foam::BlockILUCpPrecon<Type>::calcFactorization
 }
 
 
+template<class Type>
+void Foam::BlockILUCpPrecon<Type>::calcFactorization() const
+{
+    // Get upper and lower matrix factors
+    CoeffField<Type>& Lower = extBlockMatrix_.extendedLower();
+    CoeffField<Type>& Upper = extBlockMatrix_.extendedUpper();
+
+    // Calculate factorization
+    // Note: lower, diag and upper must have same type as required by the
+    // algorithm. This is handled by lowest possible promotion
+    if (preconDiag_.activeType() == blockCoeffBase::SCALAR)
+    {
+        if (Upper.activeType() == blockCoeffBase::SCALAR)
+        {
+            calcActiveTypeFactorization
+            (
+                preconDiag_.asScalar(),
+                Upper.asScalar(),
+                Lower.asScalar()
+            );
+        }
+        else if (Upper.activeType() == blockCoeffBase::LINEAR)
+        {
+            calcActiveTypeFactorization
+            (
+                preconDiag_.asLinear(), // Promotes to linear
+                Upper.asLinear(),
+                Lower.asLinear()
+            );
+        }
+        else if (Upper.activeType() == blockCoeffBase::SQUARE)
+        {
+            calcActiveTypeFactorization
+            (
+                preconDiag_.asSquare(), // Promotes to square
+                Upper.asSquare(),
+                Lower.asSquare()
+            );
+        }
+    }
+    else if (preconDiag_.activeType() == blockCoeffBase::LINEAR)
+    {
+        if (Upper.activeType() == blockCoeffBase::SCALAR)
+        {
+            calcActiveTypeFactorization
+            (
+                preconDiag_.asLinear(),
+                Upper.asLinear(), // Promotes to linear
+                Lower.asLinear()  // Promotes to linear
+            );
+        }
+        else if (Upper.activeType() == blockCoeffBase::LINEAR)
+        {
+            calcActiveTypeFactorization
+            (
+                preconDiag_.asLinear(),
+                Upper.asLinear(),
+                Lower.asLinear()
+            );
+        }
+        else if (Upper.activeType() == blockCoeffBase::SQUARE)
+        {
+            calcActiveTypeFactorization
+            (
+                preconDiag_.asSquare(), // Promotes to square
+                Upper.asSquare(),
+                Lower.asSquare()
+            );
+        }
+    }
+    else if (preconDiag_.activeType() == blockCoeffBase::SQUARE)
+    {
+        if (Upper.activeType() == blockCoeffBase::SCALAR)
+        {
+            calcActiveTypeFactorization
+            (
+                preconDiag_.asSquare(),
+                Upper.asSquare(), // Promotes to square
+                Lower.asSquare()  // Promotes to square
+            );
+        }
+        else if (Upper.activeType() == blockCoeffBase::LINEAR)
+        {
+            calcActiveTypeFactorization
+            (
+                preconDiag_.asSquare(),
+                Upper.asSquare(), // Promotes to square
+                Lower.asSquare()  // Promotes to square
+            );
+        }
+        else if (Upper.activeType() == blockCoeffBase::SQUARE)
+        {
+            calcActiveTypeFactorization
+            (
+                preconDiag_.asSquare(),
+                Upper.asSquare(),
+                Lower.asSquare()
+            );
+        }
+    }
+}
+
+
 template<class Type>
 template<class LDUType>
 void Foam::BlockILUCpPrecon<Type>::LUSubstitute
@@ -420,135 +531,9 @@ Foam::BlockILUCpPrecon<Type>::BlockILUCpPrecon
         (
             polyMesh::defaultRegion
         )
-    ),
+    )
-    zDiag_(1),
-    z_(preconDiag_.size()),
-    w_(preconDiag_.size())
 {
-    // Get upper and lower matrix factors
+    calcFactorization();
-    CoeffField<Type>& Lower = extBlockMatrix_.extendedLower();
-    CoeffField<Type>& Upper = extBlockMatrix_.extendedUpper();
-
-    // Calculate factorization
-    // Note: lower, diag and upper must have same type as required by the
-    // algorithm. This is handled by lowest possible promotion
-    if (preconDiag_.activeType() == blockCoeffBase::SCALAR)
-    {
-        if (Upper.activeType() == blockCoeffBase::SCALAR)
-        {
-            calcFactorization
-            (
-                preconDiag_.asScalar(),
-                Upper.asScalar(),
-                Lower.asScalar(),
-                zDiag_.asScalar(),
-                z_.asScalar(),
-                w_.asScalar()
-            );
-        }
-        else if (Upper.activeType() == blockCoeffBase::LINEAR)
-        {
-            calcFactorization
-            (
-                preconDiag_.asLinear(), // Promotes to linear
-                Upper.asLinear(),
-                Lower.asLinear(),
-                zDiag_.asLinear(),
-                z_.asLinear(),
-                w_.asLinear()
-            );
-        }
-        else if (Upper.activeType() == blockCoeffBase::SQUARE)
-        {
-            calcFactorization
-            (
-                preconDiag_.asSquare(), // Promotes to square
-                Upper.asSquare(),
-                Lower.asSquare(),
-                zDiag_.asSquare(),
-                z_.asSquare(),
-                w_.asSquare()
-            );
-        }
-    }
-    else if (preconDiag_.activeType() == blockCoeffBase::LINEAR)
-    {
-        if (Upper.activeType() == blockCoeffBase::SCALAR)
-        {
-            calcFactorization
-            (
-                preconDiag_.asLinear(),
-                Upper.asLinear(), // Promotes to linear
-                Lower.asLinear(), // Promotes to linear
-                zDiag_.asLinear(), // Promotes to linear
-                z_.asLinear(), // Promotes to linear
-                w_.asLinear() // Promotes to linear
-            );
-        }
-        else if (Upper.activeType() == blockCoeffBase::LINEAR)
-        {
-            calcFactorization
-            (
-                preconDiag_.asLinear(),
-                Upper.asLinear(),
-                Lower.asLinear(),
-                zDiag_.asLinear(),
-                z_.asLinear(),
-                w_.asLinear()
-            );
-        }
-        else if (Upper.activeType() == blockCoeffBase::SQUARE)
-        {
-            calcFactorization
-            (
-                preconDiag_.asSquare(), // Promotes to square
-                Upper.asSquare(),
-                Lower.asSquare(),
-                zDiag_.asSquare(),
-                z_.asSquare(),
-                w_.asSquare()
-            );
-        }
-    }
-    else if (preconDiag_.activeType() == blockCoeffBase::SQUARE)
-    {
-        if (Upper.activeType() == blockCoeffBase::SCALAR)
-        {
-            calcFactorization
-            (
-                preconDiag_.asSquare(),
-                Upper.asSquare(), // Promotes to square
-                Lower.asSquare(), // Promotes to square
-                zDiag_.asSquare(), // Promotes to square
-                z_.asSquare(), // Promotes to square
-                w_.asSquare() // Promotes to square
-            );
-        }
-        else if (Upper.activeType() == blockCoeffBase::LINEAR)
-        {
-            calcFactorization
-            (
-                preconDiag_.asSquare(),
-                Upper.asSquare(), // Promotes to square
-                Lower.asSquare(), // Promotes to square
-                zDiag_.asSquare(), // Promotes to square
-                z_.asSquare(), // Promotes to square
-                w_.asSquare() // Promotes to square
-            );
-        }
-        else if (Upper.activeType() == blockCoeffBase::SQUARE)
-        {
-            calcFactorization
-            (
-                preconDiag_.asSquare(),
-                Upper.asSquare(),
-                Lower.asSquare(),
-                zDiag_.asSquare(),
-                z_.asSquare(),
-                w_.asSquare()
-            );
-        }
-    }
 }
 
 

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

@@ -63,22 +63,13 @@ class BlockILUCpPrecon
     // Private Data
 
         //- Preconditioned diagonal
-        CoeffField<Type> preconDiag_;
+        mutable CoeffField<Type> preconDiag_;
 
         //- Fill in level
         const label p_;
 
         //- Extended lduMarix
-        extendedBlockLduMatrix<Type> extBlockMatrix_;
+        mutable extendedBlockLduMatrix<Type> extBlockMatrix_;
-
-        //- Temporary working diagonal
-        CoeffField<Type> zDiag_;
-
-        //- Temporary working row field
-        CoeffField<Type> z_;
-
-        //- Temporary working column field
-        CoeffField<Type> w_;
 
 
     // Private Member Functions
@@ -89,18 +80,18 @@ class BlockILUCpPrecon
         //- Disallow default bitwise assignment
         void operator=(const BlockILUCpPrecon&);
 
-        //- Calculate LU factorization.
+        //- Calculate active type dependant factorization.
         //  Note: both lower, diag and upper have to be of same type
         template<class LDUType>
-        void calcFactorization
+        void calcActiveTypeFactorization
         (
             Field<LDUType>& preconD,
             Field<LDUType>& extUpper,
-            Field<LDUType>& extLower,
+            Field<LDUType>& extLower
-            Field<LDUType>& zDiag,
+        ) const;
-            Field<LDUType>& z,
+
-            Field<LDUType>& w
+        //- Calculate LU factorization (constructor helper)
-        );
+        void calcFactorization() const;
 
         //- Performs forward/backward LU substitution
         //  Note: both lower, diag and upper have to be of same type

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/blockILUCpPrecons.C

@@ -39,10 +39,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-// Excluded tensor, need reduced version. VV, 3/Jul/2015.
+makeBlockPrecons(blockILUCpPrecon);
-makeBlockPrecon(blockScalarPrecon, blockILUCpPreconScalar);
-makeBlockPrecon(blockVectorPrecon, blockILUCpPreconVector);
-//makeBlockPrecon(blockTensorPrecon, blockILUCpPreconTensor);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/blockILUCpPrecons.H

@@ -42,8 +42,6 @@ SourceFiles
 #include "scalarBlockILUCpPrecon.H"
 #include "tensorBlockILUCpPrecon.H"
 
-#include "BlockILUCpPrecon.H"
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/scalarBlockILUCpPrecon.C

@@ -45,18 +45,23 @@ namespace Foam
 
 template<>
 template<>
-void BlockILUCpPrecon<scalar>::calcFactorization
+void BlockILUCpPrecon<scalar>::calcActiveTypeFactorization
 (
     scalarField& preconD,
     scalarField& extUpper,
-    scalarField& extLower,
+    scalarField& extLower
-    scalarField& zDiag,
+) const
-    scalarField& z,
-    scalarField& w
-)
 {
     if (!matrix_.diagonal())
     {
+        // Get number of rows
+        const label nRows = preconD.size();
+
+        // Allocate working fields
+        scalarField z(nRows, 0.0);
+        scalarField w(nRows, 0.0);
+        scalar zDiag = 0.0;
+
         // Get necessary const access to extended ldu addressing
         const extendedLduAddressing& addr = extBlockMatrix_.extendedLduAddr();
 
@@ -82,10 +87,6 @@ void BlockILUCpPrecon<scalar>::calcFactorization
         // Get access to working fields
         scalar* __restrict__ zPtr = z.begin();
         scalar* __restrict__ wPtr = w.begin();
-        scalar& zDiagI = zDiag[0];
-
-        // Get number of rows
-        const label nRows = preconD.size();
 
         // Define start and end face ("virtual" face when extended addressing is
         // used) of this row/column.
@@ -101,7 +102,7 @@ void BlockILUCpPrecon<scalar>::calcFactorization
             fEnd = ownStartPtr[rowI + 1];
 
             // Initialize temporary working diagonal
-            zDiagI = diagPtr[rowI];
+            zDiag = diagPtr[rowI];
 
             // Initialize temporary working row field
             for (register label faceI = fStart; faceI < fEnd; ++faceI)
@@ -132,7 +133,7 @@ void BlockILUCpPrecon<scalar>::calcFactorization
                 const label i = lPtr[losortCoeff];
 
                 // Update diagonal
-                zDiagI -= lowerPtr[losortCoeff]*upperPtr[losortCoeff];
+                zDiag -= lowerPtr[losortCoeff]*upperPtr[losortCoeff];
 
                 // Get end of row for cell i
                 const register label fEndRowi = ownStartPtr[i + 1];
@@ -154,7 +155,7 @@ void BlockILUCpPrecon<scalar>::calcFactorization
 
             // Update diagonal entry, inverting it for future use
             scalar& diagRowI = diagPtr[rowI];
-            diagRowI = 1.0/zDiagI;
+            diagRowI = 1.0/zDiag;
 
             // Index for updating L and U
             register label zwIndex;
@@ -171,7 +172,7 @@ void BlockILUCpPrecon<scalar>::calcFactorization
             }
 
             // Reset temporary working fields
-            zDiagI = 0;
+            zDiag = 0;
 
             // Only reset parts of the working fields that have been updated in
             // this step (for this row and column)
@@ -228,6 +229,18 @@ void BlockILUCpPrecon<scalar>::calcFactorization
 }
 
 
+template<>
+void BlockILUCpPrecon<scalar>::calcFactorization() const
+{
+    calcActiveTypeFactorization
+    (
+        preconDiag_.asScalar(),
+        extBlockMatrix_.extendedUpper().asScalar(),
+        extBlockMatrix_.extendedLower().asScalar()
+    );
+}
+
+
 template<>
 void BlockILUCpPrecon<scalar>::precondition
 (

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/scalarBlockILUCpPrecon.H

@@ -45,18 +45,21 @@ SourceFiles
 namespace Foam
 {
 
-// Calculate factorization
+// Calculate active type factorization
 template<>
 template<>
-void BlockILUCpPrecon<scalar>::calcFactorization
+void BlockILUCpPrecon<scalar>::calcActiveTypeFactorization
 (
     scalarField& preconD,
     scalarField& extUpper,
-    scalarField& extLower,
+    scalarField& extLower
-    scalarField& zDiag,
+) const;
-    scalarField& z,
+
-    scalarField& w
+
-);
+// Calculate factorization (constructor helper)
+template<>
+void BlockILUCpPrecon<scalar>::calcFactorization() const;
+
 
 // Precondition
 template<>

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/tensorBlockILUCpPrecon.C

@@ -43,24 +43,75 @@ Author
 namespace Foam
 {
 
-// Calculate factorization
 template<>
 template<>
-void BlockILUCpPrecon<tensor>::calcFactorization
+void BlockILUCpPrecon<tensor>::calcActiveTypeFactorization
 (
     tensorField& preconD,
     tensorField& extUpper,
-    tensorField& extLower,
+    tensorField& extLower
-    tensorField& zDiag,
+) const
-    tensorField& z,
-    tensorField& w
-)
 {
     // Decoupled version
     notImplemented("void Foam::BlockILUCpPrecon<tensor>::calcFactorization");
 }
 
 
+template<>
+void BlockILUCpPrecon<tensor>::calcFactorization() const
+{
+    // Get upper and lower matrix factors
+    CoeffField<tensor>& Lower = extBlockMatrix_.extendedLower();
+    CoeffField<tensor>& Upper = extBlockMatrix_.extendedUpper();
+
+    // Calculate factorization
+    // Note: lower, diag and upper must have same type as required by the
+    // algorithm. This is handled by lowest possible promotion
+    if (preconDiag_.activeType() == blockCoeffBase::SCALAR)
+    {
+        if (Upper.activeType() == blockCoeffBase::SCALAR)
+        {
+            calcActiveTypeFactorization
+            (
+                preconDiag_.asScalar(),
+                Upper.asScalar(),
+                Lower.asScalar()
+            );
+        }
+        else if (Upper.activeType() == blockCoeffBase::LINEAR)
+        {
+            calcActiveTypeFactorization
+            (
+                preconDiag_.asLinear(), // Promotes to linear
+                Upper.asLinear(),
+                Lower.asLinear()
+            );
+        }
+    }
+    else if (preconDiag_.activeType() == blockCoeffBase::LINEAR)
+    {
+        if (Upper.activeType() == blockCoeffBase::SCALAR)
+        {
+            calcActiveTypeFactorization
+            (
+                preconDiag_.asLinear(),
+                Upper.asLinear(), // Promotes to linear
+                Lower.asLinear()  // Promotes to linear
+            );
+        }
+        else if (Upper.activeType() == blockCoeffBase::LINEAR)
+        {
+            calcActiveTypeFactorization
+            (
+                preconDiag_.asLinear(),
+                Upper.asLinear(),
+                Lower.asLinear()
+            );
+        }
+    }
+}
+
+
 template<>
 void BlockILUCpPrecon<tensor>::precondition
 (

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/tensorBlockILUCpPrecon.H

@@ -45,20 +45,23 @@ SourceFiles
 namespace Foam
 {
 
-// Calculate factorization
+// Calculate active type factorization
 template<>
 template<>
-void BlockILUCpPrecon<tensor>::calcFactorization
+void BlockILUCpPrecon<tensor>::calcActiveTypeFactorization
 (
     tensorField& preconD,
     tensorField& extUpper,
-    tensorField& extLower,
+    tensorField& extLower
-    tensorField& zDiag,
+) const;
-    tensorField& z,
-    tensorField& w
-);
 
 
+// Calculate factorization (constructor helper)
+template<>
+void BlockILUCpPrecon<tensor>::calcFactorization() const;
+
+
+// Precondition
 template<>
 void BlockILUCpPrecon<tensor>::precondition
 (
@@ -67,6 +70,7 @@ void BlockILUCpPrecon<tensor>::precondition
 ) const;
 
 
+// Precondition transpose
 template<>
 void BlockILUCpPrecon<tensor>::preconditionT
 (

src/foam/primitives/BlockCoeff/BlockCoeff.H

@@ -126,6 +126,37 @@ public:
             }
 
 
+        // Coefficient times coefficient multiplication. Needed for BlockILUCp
+        // preconditioner. VV, 12/Jul/2015.
+
+            scalarType activeTypeMultiply
+            (
+                const scalarType& a,
+                const scalarType& b
+            ) const
+            {
+                return a*b;
+            }
+
+            linearType activeTypeMultiply
+            (
+                const linearType& a,
+                const linearType& b
+            ) const
+            {
+                return cmptMultiply(a, b);
+            }
+
+            squareType activeTypeMultiply
+            (
+                const squareType& a,
+                const squareType& b
+            ) const
+            {
+                return (a & b);
+            }
+
+
         // Transpose functions
 
             scalarType transpose(const scalarType& c) const

src/foam/primitives/BlockCoeff/DecoupledBlockCoeff.H

@@ -113,6 +113,28 @@ public:
         }
 
 
+        // Coefficient times coefficient multiplication. Needed for BlockILUCp
+        // preconditioner. VV, 12/Jul/2015.
+
+            scalarType activeTypeMultiply
+            (
+                const scalarType& a,
+                const scalarType& b
+            ) const
+            {
+                return a*b;
+            }
+
+            linearType activeTypeMultiply
+            (
+                const linearType& a,
+                const linearType& b
+            ) const
+            {
+                return cmptMultiply(a, b);
+            }
+
+
         // Inverse functions
 
             scalarType inverse(const scalarType& c) const