Minor extendedLdu matrices consistency rewrite compared to ldu matrices

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

@@ -27,12 +27,18 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+template<class Type>
+void Foam::extendedBlockLduMatrix<Type>::clearOut()
+{
+    deleteDemandDrivenData(extendedLowerPtr_);
+    deleteDemandDrivenData(extendedUpperPtr_);
+}
+
+
 template<class Type>
 void Foam::extendedBlockLduMatrix<Type>::mapOffDiagCoeffs
 (
-    const BlockLduMatrix<Type>& blockLdum,
+    const BlockLduMatrix<Type>& blockLdum
-    const label extensionLevel,
-    const polyMesh& polyMesh
 )
 {
     if (blockLdum.diagonal())
@@ -52,17 +58,22 @@ void Foam::extendedBlockLduMatrix<Type>::mapOffDiagCoeffs
         // Get reference to faceMap in extended addressing
         const unallocLabelList& faceMap = extLduAddr_.faceMap();
 
-        // Avoid assuming it's upper if the matrix is symmetric
+        // Matrix is considered symmetric if the upper is allocated and lower
-        if (blockLdum.hasUpper())
+        // is not allocated. Allocating extended upper only.
-        {
-            // 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();
@@ -71,18 +82,13 @@ void Foam::extendedBlockLduMatrix<Type>::mapOffDiagCoeffs
             // positions
             forAll (upper, faceI)
             {
-                extUpper[faceMap[faceI]] = upper[faceI];
+                activeExtUpper[faceMap[faceI]] = activeUpper[faceI];
             }
         }
         else if (upper.activeType() == blockCoeffBase::LINEAR)
         {
-            // Allocate extended lower only
+            // Helper type definition
-            extendedLowerPtr_ = new TypeCoeffField
+            typedef typename CoeffField<Type>::linearTypeField activeType;
-            (
-                extLduAddr_.extendedLowerAddr().size(),
-                pTraits<Type>::zero
-            );
-            TypeCoeffField& extLower = *extendedLowerPtr_;
 
             // Get references to fields
             const activeType& activeUpper = upper.asLinear();
@@ -90,9 +96,9 @@ void Foam::extendedBlockLduMatrix<Type>::mapOffDiagCoeffs
 
             // Copy non-zero coeffs from basic lduMatrix into corresponding
             // positions
-            forAll (lower, faceI)
+            forAll (upper, faceI)
             {
-                extLower[faceMap[faceI]] = lower[faceI];
+                activeExtUpper[faceMap[faceI]] = activeUpper[faceI];
             }
         }
         else if (upper.activeType() == blockCoeffBase::SQUARE)
@@ -146,8 +152,72 @@ void Foam::extendedBlockLduMatrix<Type>::mapOffDiagCoeffs
         // Assuming lower and upper have the same active type
         if (upper.activeType() == blockCoeffBase::SCALAR)
         {
-            extUpper[faceMap[faceI]] = upper[faceI];
+            // Helper type definition
-            extLower[faceMap[faceI]] = lower[faceI];
+            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)
+            {
+                activeExtUpper[faceMap[faceI]] = activeUpper[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
+            (
+                "void extendedBlockLduMatrix<Type>::mapOffDiagCoeffs\n"
+                "(\n"
+                "    const BlockLduMatrix<Type>& blockLdum\n"
+                ")"
+            )   << "Problem between ordinary block matrix and extended"
+                << " block matrix upper/lower coeffs type morphing."
+                << abort(FatalError);
         }
     }
 }
@@ -191,15 +261,7 @@ Foam::extendedBlockLduMatrix<Type>::extendedBlockLduMatrix
 template<class Type>
 Foam::extendedBlockLduMatrix<Type>::~extendedBlockLduMatrix()
 {
-    if (extendedLowerPtr_)
+    clearOut();
-    {
-        delete extendedLowerPtr_;
-    }
-
-    if (extendedUpperPtr_)
-    {
-        delete extendedUpperPtr_;
-    }
 }
 
 

src/lduSolvers/lduPrecon/ILUCp/ILUCp.C

@@ -96,8 +96,8 @@ void Foam::ILUCp::calcFactorization()
         // Get number of rows
         const label nRows = preconDiag_.size();
 
-        // Define start and end face ("virtual" face when extended addressing is
+        // Define start and end face ("virtual" face when extended addressing
-        // used) of this row/column, and number of non zero off diagonal entries
+        // is used) of this row/column.
         register label fStart, fEnd, fLsrStart, fLsrEnd;
 
         // Crout LU factorization
@@ -244,11 +244,10 @@ Foam::ILUCp::ILUCp
     extMatrix_
     (
         matrix,
-        p_,
+        matrix.mesh().lduAddr().extendedAddr
-        matrix.mesh().thisDb().parent().lookupObject
+        (
-        <
+            readLabel(dict.lookup("fillInLevel"))
-            polyMesh
+        )
-        >(polyMesh::defaultRegion)
     ),
     preconDiag_(matrix_.diag()),
     zDiag_(0),