Second version of BlockILUCp preconditioner

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

@@ -30,7 +30,7 @@ License
 template<class Type>
 Foam::extendedBlockLduMatrix<Type>::extendedBlockLduMatrix
 (
-    const blockLduMatrix<Type>& blockLdum,
+    const BlockLduMatrix<Type>& blockLdum,
     const label extensionLevel,
     const polyMesh& polyMesh
 )
@@ -68,24 +68,74 @@ Foam::extendedBlockLduMatrix<Type>::extendedBlockLduMatrix
         const unallocLabelList& faceMap = extLduAddr_.faceMap();
 
         // Avoid assuming it's upper if the matrix is symmetric
-        if (blockLdum.hasUpper())
+        if (blockLdum.thereIsUpper())
         {
             // Allocate extended upper only
             extendedUpperPtr_ = new TypeCoeffField
             (
-                extLduAddr_.extendedUpperAddr().size(),
+                extLduAddr_.extendedUpperAddr().size()
-                pTraits<Type>::zero
             );
             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<Type>::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)
                 {
-                extUpper[faceMap[faceI]] = upper[faceI];
+                    activeExtUpper[faceMap[faceI]] = activeUpper[faceI];
+                }
+            }
+            else if (upper.activeType() == blockCoeffBase::LINEAR)
+            {
+                // Helper type definition
+                typedef typename CoeffField<Type>::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 if (upper.activeType() == blockCoeffBase::SQUARE)
+            {
+                // Helper type definition
+                typedef typename CoeffField<Type>::squareTypeField activeType;
+
+                // Get references to fields
+                const activeType& activeUpper = upper.asSquare();
+                activeType& activeExtUpper = extUpper.asSquare();
+
+                // 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
@@ -93,19 +143,69 @@ Foam::extendedBlockLduMatrix<Type>::extendedBlockLduMatrix
             // Allocate extended lower only
             extendedLowerPtr_ = new TypeCoeffField
             (
-                extLduAddr_.extendedLowerAddr().size(),
+                extLduAddr_.extendedLowerAddr().size()
-                pTraits<Type>::zero
             );
             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<Type>::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)
                 {
-                extLower[faceMap[faceI]] = lower[faceI];
+                    activeExtLower[faceMap[faceI]] = activeLower[faceI];
+                }
+            }
+            else if (lower.activeType() == blockCoeffBase::LINEAR)
+            {
+                // Helper type definition
+                typedef typename CoeffField<Type>::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 if (lower.activeType() == blockCoeffBase::SQUARE)
+            {
+                // Helper type definition
+                typedef typename CoeffField<Type>::squareTypeField activeType;
+
+                // Get references to fields
+                const activeType& activeLower = lower.asSquare();
+                activeType& activeExtLower = extLower.asSquare();
+
+                // 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);
             }
         }
     }
@@ -118,22 +218,82 @@ Foam::extendedBlockLduMatrix<Type>::extendedBlockLduMatrix
         const label nExtFaces = extLduAddr_.extendedUpperAddr().size();
 
         // Allocate extended upper and lower
-        extendedUpperPtr_ = new TypeCoeffField(nExtFaces, pTraits<Type>::zero);
+        extendedUpperPtr_ = new TypeCoeffField(nExtFaces);
         TypeCoeffField& extUpper = *extendedUpperPtr_;
 
-        extendedLowerPtr_ = new TypeCoeffField(nExtFaces, pTraits<Type>::zero);
+        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<Type>::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)
             {
-            extUpper[faceMap[faceI]] = upper[faceI];
+                activeExtUpper[faceMap[faceI]] = activeUpper[faceI];
-            extLower[faceMap[faceI]] = lower[faceI];
+                activeExtLower[faceMap[faceI]] = activeLower[faceI];
+            }
+        }
+        else if (upper.activeType() == blockCoeffBase::LINEAR)
+        {
+            // Helper type definition
+            typedef typename CoeffField<Type>::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 if (upper.activeType() == blockCoeffBase::SQUARE)
+        {
+            // Helper type definition
+            typedef typename CoeffField<Type>::squareTypeField activeType;
+
+            // Get references to fields
+            const activeType& activeUpper = upper.asSquare();
+            activeType& activeExtUpper = extUpper.asSquare();
+            const activeType& activeLower = lower.asSquare();
+            activeType& activeExtLower = extLower.asSquare();
+
+            // 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);
         }
     }
 }
@@ -172,8 +332,7 @@ Foam::extendedBlockLduMatrix<Type>::extendedLower()
         {
             extendedLowerPtr_ = new TypeCoeffField
             (
-                extLduAddr_.extendedLowerAddr().size(),
+                extLduAddr_.extendedLowerAddr().size()
-                pTraits<Type>::zero
             );
         }
     }
@@ -196,8 +355,7 @@ Foam::extendedBlockLduMatrix<Type>::extendedUpper()
         {
             extendedUpperPtr_ = new TypeCoeffField
             (
-                extLduAddr_.extendedUpperAddr().size(),
+                extLduAddr_.extendedUpperAddr().size()
-                pTraits<Type>::zero
             );
         }
     }

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

@@ -58,7 +58,7 @@ class extendedBlockLduMatrix
 {
     // Public data type
 
-        typedef typename BlockLduMatrix<Type>::TypeCoeffField TypeCoeffField;
+        typedef CoeffField<Type> TypeCoeffField;
 
 
 private:
@@ -145,6 +145,12 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef NoRepository
+#   include "extendedBlockLduMatrix.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //

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

@@ -54,7 +54,7 @@ void Foam::BlockILUCpPrecon<Type>::calcFactorization
     if (!this->matrix_.diagonal())
     {
         // Create multiplication function object
-        typename BlockCoeff<Type>::multiply mult;
+        typename BlockCoeff<LDUType>::multiply mult;
 
         // Get necessary const access to extended ldu addressing
         const extendedLduAddressing& addr = extBlockMatrix_.extendedLduAddr();
@@ -98,7 +98,6 @@ 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];
 
@@ -172,7 +171,7 @@ void Foam::BlockILUCpPrecon<Type>::calcFactorization
 
             // Update diagonal entry, inverting it for future use
             LDUType& diagRowI = diagPtr[rowI];
-            diagRowI = mult.inv(zDiagI);
+            diagRowI = mult.inverse(zDiagI);
 
             // Index for updating L and U
             register label zwIndex;

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

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

src/foam/matrices/blockLduMatrix/BlockLduSolvers/blockVectorNSolvers.C

@@ -31,6 +31,7 @@ License
 #include "blockDiagonalPrecons.H"
 #include "blockGaussSeidelPrecons.H"
 #include "BlockCholeskyPrecon.H"
+#include "BlockILUCpPrecon.H"
 
 #include "blockLduSmoothers.H"
 #include "blockGaussSeidelSmoothers.H"
@@ -83,6 +84,9 @@ makeBlockPrecon(block##Type##Precon, block##Type##GaussSeidelPrecon);         \
 typedef BlockCholeskyPrecon<type > block##Type##CholeskyPrecon;               \
 makeBlockPrecon(block##Type##Precon, block##Type##CholeskyPrecon);            \
                                                                               \
+typedef BlockILUCpPrecon<type > block##Type##ILUCpPrecon;                     \
+makeBlockPrecon(block##Type##Precon, block##Type##ILUCpPrecon);               \
+                                                                              \
 /* Smoothers */                                                               \
 typedef BlockLduSmoother<type > block##Type##Smoother;                        \
 defineNamedTemplateTypeNameAndDebug(block##Type##Smoother, 0);                \