From 86fc3380465d7f8158892773fcb390b1861f00d5 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Wed, 28 Oct 2015 10:08:07 +0000
Subject: [PATCH] Feature ILUCp and BlockILUCp preconditioning, extended addr
 update and clean-up

---
 .../expandContract/ExpandTensorNField.C       |   6 +-
 .../expandContract/ExpandTensorNField.H       |   5 +-
 .../fields/expandContract/expandTensorField.C |  19 ++
 .../fields/expandContract/expandTensorField.H |   3 +
 .../BlockAmg/coarseBlockAmgLevel.C            |   9 +-
 .../BlockAmg/fineBlockAmgLevel.C              |   9 +-
 .../BlockLduMatrix/BlockLduMatrix.H           |   5 +-
 .../BlockLduMatrix/BlockLduMatrixOperations.C | 207 +++++++++++++-----
 .../extendedBlockLduMatrix.C                  |  42 ++--
 .../extendedBlockLduMatrix.H                  |   6 +-
 .../sphericalTensorExtendedBlockLduMatrix.C   |  27 ++-
 .../symmTensorExtendedBlockLduMatrix.C        |  32 ++-
 .../tensorExtendedBlockLduMatrix.C            |  36 ++-
 .../BlockGaussSeidelPrecon.C                  |   3 +-
 .../BlockILUCpPrecon/BlockILUCpPrecon.C       |  20 +-
 .../BlockILUCpPrecon/BlockILUCpPrecon.H       |   9 +-
 .../extendedLduAddressing.C                   |  79 ++++---
 .../extendedLduAddressing.H                   |  50 ++---
 .../lduMatrix/lduAddressing/lduAddressing.C   |  36 +++
 .../lduMatrix/lduAddressing/lduAddressing.H   |  24 +-
 .../extendedLduMatrix/extendedLduMatrix.C     |  21 +-
 .../extendedLduMatrix/extendedLduMatrix.H     |  14 +-
 src/foam/primitives/VectorN/TensorNI.H        |  42 +++-
 src/foam/primitives/VectorN/VectorN.H         |   2 +
 src/foam/primitives/VectorN/VectorNI.H        |   6 +
 .../VectorN/expandContract/ExpandTensorN.H    |   3 +
 .../VectorN/expandContract/ExpandTensorNI.H   |  86 +++++++-
 .../primitives/expandContract/expandTensor.H  |  48 ++++
 src/lduSolvers/lduPrecon/ILUCp/ILUCp.C        |  26 +--
 src/lduSolvers/lduPrecon/ILUCp/ILUCp.H        |  11 +-
 30 files changed, 632 insertions(+), 254 deletions(-)

diff --git a/src/foam/fields/Fields/VectorNFields/expandContract/ExpandTensorNField.C b/src/foam/fields/Fields/VectorNFields/expandContract/ExpandTensorNField.C
index 51ba1d848..248a87247 100644
--- a/src/foam/fields/Fields/VectorNFields/expandContract/ExpandTensorNField.C
+++ b/src/foam/fields/Fields/VectorNFields/expandContract/ExpandTensorNField.C
@@ -77,7 +77,11 @@ UNARY_FUNCTION(sphericalTensorType, cmptType, expandScalar)                   \
                                                                               \
 UNARY_FUNCTION(tensorType, vectorType, expandLinear)                          \
 UNARY_FUNCTION(diagTensorType, vectorType, expandLinear)                      \
-UNARY_FUNCTION(sphericalTensorType, vectorType, expandLinear)
+UNARY_FUNCTION(sphericalTensorType, vectorType, expandLinear)                 \
+                                                                              \
+UNARY_FUNCTION(vectorType, tensorType, sumToDiag)                             \
+UNARY_FUNCTION(vectorType, tensorType, sumMagToDiag)
+
 
 namespace Foam
 {
diff --git a/src/foam/fields/Fields/VectorNFields/expandContract/ExpandTensorNField.H b/src/foam/fields/Fields/VectorNFields/expandContract/ExpandTensorNField.H
index 2c4b7da2e..71b123cae 100644
--- a/src/foam/fields/Fields/VectorNFields/expandContract/ExpandTensorNField.H
+++ b/src/foam/fields/Fields/VectorNFields/expandContract/ExpandTensorNField.H
@@ -65,7 +65,10 @@ UNARY_FUNCTION(sphericalTensorType, cmptType, expandScalar)                   \
                                                                               \
 UNARY_FUNCTION(tensorType, vectorType, expandLinear)                          \
 UNARY_FUNCTION(diagTensorType, vectorType, expandLinear)                      \
-UNARY_FUNCTION(sphericalTensorType, vectorType, expandLinear)
+UNARY_FUNCTION(sphericalTensorType, vectorType, expandLinear)                 \
+                                                                              \
+UNARY_FUNCTION(vectorType, tensorType, sumToDiag)                             \
+UNARY_FUNCTION(vectorType, tensorType, sumMagToDiag)
 
 
 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
diff --git a/src/foam/fields/expandContract/expandTensorField.C b/src/foam/fields/expandContract/expandTensorField.C
index 7f4d9f946..fb0698104 100644
--- a/src/foam/fields/expandContract/expandTensorField.C
+++ b/src/foam/fields/expandContract/expandTensorField.C
@@ -85,6 +85,25 @@ void expandLinear(Field<tensor>& res, const UList<vector>& f)
     }
 }
 
+
+void sumToDiag(Field<vector>& res, const UList<tensor>& f)
+{
+    forAll (res, i)
+    {
+        sumToDiag(res[i], f[i]);
+    }
+}
+
+
+void sumMagToDiag(Field<vector>& res, const UList<tensor>& f)
+{
+    forAll (res, i)
+    {
+        sumMagToDiag(res[i], f[i]);
+    }
+}
+
+
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/foam/fields/expandContract/expandTensorField.H b/src/foam/fields/expandContract/expandTensorField.H
index 1ae41828f..2cfc568b8 100644
--- a/src/foam/fields/expandContract/expandTensorField.H
+++ b/src/foam/fields/expandContract/expandTensorField.H
@@ -54,6 +54,9 @@ void expandScalar(Field<vector>& res, const UList<scalar>& f);
 void expandScalar(Field<tensor>& res, const UList<scalar>& f);
 void expandLinear(Field<tensor>& res, const UList<vector>& f);
 
+void sumToDiag(Field<vector>& res, const UList<tensor>& f);
+void sumMagToDiag(Field<vector>& res, const UList<tensor>& f);
+
 } // End namespace Foam
 
 #endif
diff --git a/src/foam/matrices/blockLduMatrix/BlockAmg/coarseBlockAmgLevel.C b/src/foam/matrices/blockLduMatrix/BlockAmg/coarseBlockAmgLevel.C
index 2f7cd9d72..5ea4a85c8 100644
--- a/src/foam/matrices/blockLduMatrix/BlockAmg/coarseBlockAmgLevel.C
+++ b/src/foam/matrices/blockLduMatrix/BlockAmg/coarseBlockAmgLevel.C
@@ -76,12 +76,15 @@ Foam::coarseBlockAmgLevel<Type>::coarseBlockAmgLevel
         BlockLduSmoother<Type>::New
         (
             matrixPtr_,
-            dict,
-            "coarseSmoother"
+            dict
+//             ,"coarseSmoother"
         )
     ),
     Ax_()
-{}
+{
+    Info<< "Coarse AMG level check" << endl;
+    matrixPtr_->check();
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
diff --git a/src/foam/matrices/blockLduMatrix/BlockAmg/fineBlockAmgLevel.C b/src/foam/matrices/blockLduMatrix/BlockAmg/fineBlockAmgLevel.C
index 14bfbc773..691649a37 100644
--- a/src/foam/matrices/blockLduMatrix/BlockAmg/fineBlockAmgLevel.C
+++ b/src/foam/matrices/blockLduMatrix/BlockAmg/fineBlockAmgLevel.C
@@ -71,12 +71,15 @@ Foam::fineBlockAmgLevel<Type>::fineBlockAmgLevel
         BlockLduSmoother<Type>::New
         (
             matrix,
-            dict,
-            "fineSmoother"
+            dict
+//             ,"fineSmoother"
         )
     ),
     Ax_()
-{}
+{
+    Info<< "Fine AMG level check" << endl;
+    matrix_.check();
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrix.H b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrix.H
index ec91d6e89..9a1722593 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrix.H
+++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrix.H
@@ -227,9 +227,8 @@ public:
         BlockLduMatrix(BlockLduMatrix<Type>&, bool reUse);
 
 
-    // Destructor
-
-        virtual ~BlockLduMatrix();
+    //- Destructor
+    virtual ~BlockLduMatrix();
 
 
     // Member functions
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixOperations.C b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixOperations.C
index dcdd3b3df..be8033239 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixOperations.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixOperations.C
@@ -307,9 +307,14 @@ void Foam::BlockLduMatrix<Type>::check() const
     typedef typename TypeCoeffField::linearTypeField linearTypeField;
     typedef typename TypeCoeffField::squareTypeField squareTypeField;
 
-    // Copy the diagonal
-    TypeCoeffField DiagCopy(this->diag().size());
-//     TypeCoeffField DiagCopy = this->diag();
+    // TODO: Account for coupled boundary patch coeffs
+    // HJ, 29/Oct/2015
+
+    // Get original diag
+    const TypeCoeffField& d = this->diag();
+
+    // Copy the diagonal.  It is initialised to zero
+    TypeCoeffField SumOffDiag(d.size());
 
     const unallocLabelList& l = lduAddr().lowerAddr();
     const unallocLabelList& u = lduAddr().upperAddr();
@@ -323,67 +328,110 @@ void Foam::BlockLduMatrix<Type>::check() const
         if
         (
             Upper.activeType() == blockCoeffBase::SQUARE
-         || DiagCopy.activeType() == blockCoeffBase::SQUARE
+         || d.activeType() == blockCoeffBase::SQUARE
         )
         {
+            // Note: For a square coefficient, the matrix needs to be analysed
+            // row-by-row, including the contribution of the off-diagonal
+            // elements in the diagonal matrix
+
+            // Get result as linear
+            linearTypeField& activeSumOffDiag = SumOffDiag.asLinear();
+
+            // Do diagonal elements
+
+            // Create the diagonal matrix element without the diagonal
+            const squareTypeField& activeDiag = d.asSquare();
+
+            linearTypeField diagDiag(activeDiag.size());
+            squareTypeField luDiag(activeDiag.size());
+
+            // Take out the diagonal from the diag as a linear type
+            contractLinear(diagDiag, activeDiag);
+
+            // Expand diagonal only to full square type and store into luDiag
+            expandLinear(luDiag, diagDiag);
+
+            // Keep only off-diagonal in luDiag
+            luDiag = activeDiag - luDiag;
+
+            sumMagToDiag(activeSumOffDiag, luDiag);
+
+            // Do the off-diagonal elements by collapsing each row
+            // into the linear form
+
             const squareTypeField& activeUpper = Upper.asSquare();
-            squareTypeField& activeDiagCopy = DiagCopy.asSquare();
 
             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
             {
-                activeDiagCopy[l[coeffI]] += activeUpper[coeffI].T();
-                activeDiagCopy[u[coeffI]] += activeUpper[coeffI];
+                activeSumOffDiag[l[coeffI]] +=
+                    sumMagToDiag(activeUpper[coeffI].T());
+
+                activeSumOffDiag[u[coeffI]] +=
+                    sumMagToDiag(activeUpper[coeffI]);
             }
 
-            Info<< "void BlockLduMatrix<Type>::check() const : "
-                << "Symmetric matrix: raw matrix difference: "
-                << sum(mag(activeDiagCopy))
-                << " scaled: "
-                << sum(mag(activeDiagCopy))/sum(mag(this->diag().asSquare()))
+            // Check diagonal dominance
+
+            diagDiag = cmptMag(diagDiag);
+
+            // Divide diagonal with sum of off-diagonals
+            cmptDivide(diagDiag, diagDiag, activeSumOffDiag);
+
+            InfoIn("void BlockLduMatrix<Type>::check() const)")
+                << "Symmetric matrix: " << activeDiag.size()
+                << " diagonal dominance sym square: "
+                << Foam::min(diagDiag)
                 << endl;
         }
         else if
         (
             Upper.activeType() == blockCoeffBase::LINEAR
-         || DiagCopy.activeType() == blockCoeffBase::LINEAR
+         || d.activeType() == blockCoeffBase::LINEAR
         )
         {
+            const linearTypeField& activeDiag = d.asLinear();
+
             const linearTypeField& activeUpper = Upper.asLinear();
-            linearTypeField& activeDiagCopy = DiagCopy.asLinear();
+            linearTypeField& activeSumOffDiag = SumOffDiag.asLinear();
 
             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
             {
-                activeDiagCopy[l[coeffI]] += activeUpper[coeffI];
-                activeDiagCopy[u[coeffI]] += activeUpper[coeffI];
+                activeSumOffDiag[l[coeffI]] += cmptMag(activeUpper[coeffI]);
+                activeSumOffDiag[u[coeffI]] += cmptMag(activeUpper[coeffI]);
             }
 
-            Info<< "void BlockLduMatrix<Type>::check() const : "
-                << "Symmetric matrix: raw matrix difference: "
-                << sum(mag(activeDiagCopy))
-                << " scaled: "
-                << sum(mag(activeDiagCopy))/sum(mag(this->diag().asLinear()))
+            linearTypeField diagDiag = cmptMag(activeDiag);
+
+            cmptDivide(diagDiag, diagDiag, activeSumOffDiag);
+
+            InfoIn("void BlockLduMatrix<Type>::check() const)")
+                << "Symmetric matrix: " << activeDiag.size()
+                << " diagonal dominance sym linear: "
+                << Foam::min(diagDiag)
                 << endl;
         }
         else if
         (
             Upper.activeType() == blockCoeffBase::SCALAR
-         || DiagCopy.activeType() == blockCoeffBase::SCALAR
+         || d.activeType() == blockCoeffBase::SCALAR
         )
         {
+            const scalarTypeField& activeDiag = d.asScalar();
+
             const scalarTypeField& activeUpper = Upper.asScalar();
-            scalarTypeField& activeDiagCopy = DiagCopy.asScalar();
+            scalarTypeField& activeSumOffDiag = SumOffDiag.asScalar();
 
             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
             {
-                activeDiagCopy[l[coeffI]] += activeUpper[coeffI];
-                activeDiagCopy[u[coeffI]] += activeUpper[coeffI];
+                activeSumOffDiag[l[coeffI]] += cmptMag(activeUpper[coeffI]);
+                activeSumOffDiag[u[coeffI]] += cmptMag(activeUpper[coeffI]);
             }
 
-            Info<< "void BlockLduMatrix<Type>::check() const : "
-                << "Symmetric matrix: raw matrix difference: "
-                << sum(mag(activeDiagCopy))
-                << " scaled: "
-                << sum(mag(activeDiagCopy))/sum(mag(this->diag().asScalar()))
+            InfoIn("void BlockLduMatrix<Type>::check() const)")
+                << "Symmetric matrix: " << activeDiag.size()
+                << " diagonal dominance sym scalar: "
+                << Foam::min(mag(activeDiag)/activeSumOffDiag)
                 << endl;
         }
     }
@@ -398,78 +446,121 @@ void Foam::BlockLduMatrix<Type>::check() const
         (
             Lower.activeType() == blockCoeffBase::SQUARE
          || Upper.activeType() == blockCoeffBase::SQUARE
-         || DiagCopy.activeType() == blockCoeffBase::SQUARE
+         || d.activeType() == blockCoeffBase::SQUARE
         )
         {
+            // Note: For a square coefficient, the matrix needs to be analysed
+            // row-by-row, including the contribution of the off-diagonal
+            // elements in the diagonal matrix
+
+            // Get result as linear
+            linearTypeField& activeSumOffDiag = SumOffDiag.asLinear();
+
+            // Do diagonal elements
+
+            // Create the diagonal matrix element without the diagonal
+            const squareTypeField& activeDiag = d.asSquare();
+
+            linearTypeField diagDiag(activeDiag.size());
+            squareTypeField luDiag(activeDiag.size());
+
+            // Take out the diagonal from the diag as a linear type
+            contractLinear(diagDiag, activeDiag);
+
+            // Expand diagonal only to full square type and store into luDiag
+            expandLinear(luDiag, diagDiag);
+
+            // Keep only off-diagonal in luDiag
+            luDiag = activeDiag - luDiag;
+
+            sumMagToDiag(activeSumOffDiag, luDiag);
+
+            // Do the off-diagonal elements by collapsing each row
+            // into the linear form
+
             const squareTypeField& activeLower = Lower.asSquare();
             const squareTypeField& activeUpper = Upper.asSquare();
-            squareTypeField& activeDiagCopy = DiagCopy.asSquare();
 
             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
             {
-                activeDiagCopy[l[coeffI]] += activeLower[coeffI];
-                activeDiagCopy[u[coeffI]] += activeUpper[coeffI];
+                activeSumOffDiag[l[coeffI]] +=
+                    sumMagToDiag(activeLower[coeffI]);
+
+                activeSumOffDiag[u[coeffI]] +=
+                    sumMagToDiag(activeUpper[coeffI]);
             }
 
-            Info<< "void BlockLduMatrix<Type>::check() const : "
-                << "Asymmetric matrix: raw matrix difference: "
-                << sum(mag(activeDiagCopy))
-                << " scaled: "
-                << sum(mag(activeDiagCopy))/sum(mag(this->diag().asSquare()))
+            // Check diagonal dominance
+
+            diagDiag = cmptMag(diagDiag);
+
+            // Divide diagonal with sum of off-diagonals
+            cmptDivide(diagDiag, diagDiag, activeSumOffDiag);
+
+            InfoIn("void BlockLduMatrix<Type>::check() const)")
+                << "Asymmetric matrix: " << activeDiag.size()
+                << " diagonal dominance assym square: "
+                << Foam::min(diagDiag)
                 << endl;
         }
         else if
         (
             Lower.activeType() == blockCoeffBase::LINEAR
          || Upper.activeType() == blockCoeffBase::LINEAR
-         || DiagCopy.activeType() == blockCoeffBase::LINEAR
+         || d.activeType() == blockCoeffBase::LINEAR
         )
         {
+            const linearTypeField& activeDiag = d.asLinear();
+
             const linearTypeField& activeLower = Lower.asLinear();
             const linearTypeField& activeUpper = Upper.asLinear();
-            linearTypeField& activeDiagCopy = DiagCopy.asLinear();
+            linearTypeField& activeSumOffDiag = SumOffDiag.asLinear();
 
             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
             {
-                activeDiagCopy[l[coeffI]] += activeLower[coeffI];
-                activeDiagCopy[u[coeffI]] += activeUpper[coeffI];
+                activeSumOffDiag[l[coeffI]] += cmptMag(activeLower[coeffI]);
+                activeSumOffDiag[u[coeffI]] += cmptMag(activeUpper[coeffI]);
             }
 
-            Info<< "void BlockLduMatrix<Type>::check() const : "
-                << "Asymmetric matrix: raw matrix difference: "
-                << sum(mag(activeDiagCopy))
-                << " scaled: "
-                << sum(mag(activeDiagCopy))/sum(mag(this->diag().asLinear()))
+            linearTypeField diagDiag = cmptMag(activeDiag);
+
+            cmptDivide(diagDiag, diagDiag, activeSumOffDiag);
+
+            InfoIn("void BlockLduMatrix<Type>::check() const)")
+                << "Asymmetric matrix: " << activeDiag.size()
+                << " diagonal dominance assym linear: "
+                << Foam::min(diagDiag)
                 << endl;
         }
         else if
         (
             Lower.activeType() == blockCoeffBase::SCALAR
          || Upper.activeType() == blockCoeffBase::SCALAR
-         || DiagCopy.activeType() == blockCoeffBase::SCALAR
+         || d.activeType() == blockCoeffBase::SCALAR
         )
         {
+            const scalarTypeField& activeDiag = d.asScalar();
+
             const scalarTypeField& activeLower = Lower.asScalar();
             const scalarTypeField& activeUpper = Upper.asScalar();
-            scalarTypeField& activeDiagCopy = DiagCopy.asScalar();
+            scalarTypeField& activeSumOffDiag = SumOffDiag.asScalar();
 
             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
             {
-                activeDiagCopy[l[coeffI]] += activeLower[coeffI];
-                activeDiagCopy[u[coeffI]] += activeUpper[coeffI];
+                activeSumOffDiag[l[coeffI]] += cmptMag(activeLower[coeffI]);
+                activeSumOffDiag[u[coeffI]] += cmptMag(activeUpper[coeffI]);
             }
 
-            Info<< "void BlockLduMatrix<Type>::check() const : "
-                << "Asymmetric matrix: raw matrix difference: "
-                << sum(mag(activeDiagCopy))
-                << " scaled: "
-                << sum(mag(activeDiagCopy))/sum(mag(this->diag().asScalar()))
+            InfoIn("void BlockLduMatrix<Type>::check() const)")
+                << "Asymmetric matrix: "  << activeDiag.size()
+                << " diagonal dominance assym scalar: "
+                << Foam::min(mag(activeDiag)/activeSumOffDiag)
                 << endl;
         }
     }
     else
     {
-        Info<< "void BlockLduMatrix<Type>::check() const : "
+        InfoIn("void BlockLduMatrix<Type>::check() const)")
             << "Diagonal matrix" << endl;
     }
 }
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.C b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.C
index 42981b44b..ff209c3d5 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.C
@@ -35,8 +35,13 @@ void Foam::extendedBlockLduMatrix<Type>::mapOffDiagCoeffs
 {
     if (blockLdum.diagonal())
     {
-        WarningIn("extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)")
-            << "Attempted to create extended lower/upper coeffs for block "
+        WarningIn
+        (
+            "void extendedBlockLduMatrix<Type>::mapOffDiagCoeffs\n"
+            "(\n"
+            "    const BlockLduMatrix<Type>& blockLdum\n"
+            ")"
+        )   << "Attempted to create extended lower/upper coeffs for block "
             << "matrix that is diagonal."
             << nl << endl;
     }
@@ -108,7 +113,10 @@ void Foam::extendedBlockLduMatrix<Type>::mapOffDiagCoeffs
         {
             FatalErrorIn
             (
-                "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
+                "void extendedBlockLduMatrix<Type>::mapOffDiagCoeffs\n"
+                "(\n"
+                "    const BlockLduMatrix<Type>& blockLdum\n"
+                ")"
             )   << "Problem between ordinary block matrix and extended"
                 << " block matrix upper coeffs type morphing."
                 << abort(FatalError);
@@ -195,7 +203,10 @@ void Foam::extendedBlockLduMatrix<Type>::mapOffDiagCoeffs
         {
             FatalErrorIn
             (
-                "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
+                "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);
@@ -210,27 +221,24 @@ template<class Type>
 Foam::extendedBlockLduMatrix<Type>::extendedBlockLduMatrix
 (
     const BlockLduMatrix<Type>& blockLdum,
-    const label extensionLevel,
-    const polyMesh& polyMesh
+    const extendedLduAddressing& extLduAddr
 )
 :
     basicBlockLduMatrix_(blockLdum),
-    extLduAddr_
-    (
-        extendedLduAddressing::New
-        (
-            polyMesh,
-            blockLdum.lduAddr(),
-            extensionLevel
-        )
-    ),
+    extLduAddr_(extLduAddr),
     extendedLowerPtr_(NULL),
     extendedUpperPtr_(NULL)
 {
     if (debug)
     {
-        Info<< "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&) :"
-               "Constructing extendedBlockLduMatrix."
+        InfoIn
+        (
+            "extendedBlockLduMatrix<Type>::extendedBlockLduMatrix\n"
+            "(\n"
+            "    const BlockLduMatrix<Type>& blockLdum,\n"
+            "    const extendedLduAddressing& extLduAddr\n"
+            ")"
+        )   << "Constructing extendedBlockLduMatrix."
             << endl;
     }
 
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.H b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.H
index ded28183d..0b9dfa6b9 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.H
+++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/extendedBlockLduMatrix.H
@@ -39,9 +39,8 @@ SourceFiles
 #define extendedBlockLduMatrix_H
 
 #include "extendedLduAddressing.H"
-#include "polyMesh.H"
-#include "className.H"
 #include "BlockLduMatrix.H"
+#include "className.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -103,8 +102,7 @@ public:
         extendedBlockLduMatrix
         (
             const BlockLduMatrix<Type>&,
-            const label,
-            const polyMesh&
+            const extendedLduAddressing&
         );
 
 
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/sphericalTensorExtendedBlockLduMatrix.C b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/sphericalTensorExtendedBlockLduMatrix.C
index 720f8f9ec..87796e29c 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/sphericalTensorExtendedBlockLduMatrix.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/sphericalTensorExtendedBlockLduMatrix.C
@@ -41,8 +41,13 @@ void Foam::extendedBlockLduMatrix<Foam::sphericalTensor>::mapOffDiagCoeffs
 {
     if (blockLdum.diagonal())
     {
-        WarningIn("extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)")
-            << "Attempted to create extended lower/upper coeffs for block "
+        WarningIn
+        (
+            "void extendedBlockLduMatrix<sphericalTensor>::mapOffDiagCoeffs\n"
+            "(\n"
+            "    const BlockLduMatrix<sphericalTensor>& blockLdum\n"
+            ")"
+        )   << "Attempted to create extended lower/upper coeffs for block "
             << "matrix that is diagonal."
             << nl << endl;
     }
@@ -102,7 +107,11 @@ void Foam::extendedBlockLduMatrix<Foam::sphericalTensor>::mapOffDiagCoeffs
             {
                 FatalErrorIn
                 (
-                    "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
+                    "void extendedBlockLduMatrix<sphericalTensor>::"
+                    "mapOffDiagCoeffs\n"
+                    "(\n"
+                    "    const BlockLduMatrix<sphericalTensor>& blockLdum\n"
+                    ")"
                 )   << "Problem between ordinary block matrix and extended"
                     << " block matrix upper coeffs type morphing."
                     << abort(FatalError);
@@ -158,7 +167,11 @@ void Foam::extendedBlockLduMatrix<Foam::sphericalTensor>::mapOffDiagCoeffs
             {
                 FatalErrorIn
                 (
-                    "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
+                    "void extendedBlockLduMatrix<sphericalTensor>::"
+                    "mapOffDiagCoeffs\n"
+                    "(\n"
+                    "    const BlockLduMatrix<sphericalTensor>& blockLdum\n"
+                    ")"
                 )   << "Problem between ordinary block matrix and extended"
                     << " block matrix lower coeffs type morphing."
                     << abort(FatalError);
@@ -229,7 +242,11 @@ void Foam::extendedBlockLduMatrix<Foam::sphericalTensor>::mapOffDiagCoeffs
         {
             FatalErrorIn
             (
-                "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
+                "void extendedBlockLduMatrix<sphericalTensor>::"
+                "mapOffDiagCoeffs\n"
+                "(\n"
+                "    const BlockLduMatrix<sphericalTensor>& blockLdum\n"
+                ")"
             )   << "Problem between ordinary block matrix and extended"
                 << " block matrix upper/lower coeffs type morphing."
                 << abort(FatalError);
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/symmTensorExtendedBlockLduMatrix.C b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/symmTensorExtendedBlockLduMatrix.C
index 63715b2d3..75296b79d 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/symmTensorExtendedBlockLduMatrix.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/symmTensorExtendedBlockLduMatrix.C
@@ -41,8 +41,13 @@ void Foam::extendedBlockLduMatrix<Foam::symmTensor>::mapOffDiagCoeffs
 {
     if (blockLdum.diagonal())
     {
-        WarningIn("extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)")
-            << "Attempted to create extended lower/upper coeffs for block "
+        WarningIn
+        (
+            "void extendedBlockLduMatrix<symmTensor>::mapOffDiagCoeffs\n"
+            "(\n"
+            "    const BlockLduMatrix<symmTensor>& blockLdum\n"
+            ")"
+        )   << "Attempted to create extended lower/upper coeffs for block "
             << "matrix that is diagonal."
             << nl << endl;
     }
@@ -102,7 +107,11 @@ void Foam::extendedBlockLduMatrix<Foam::symmTensor>::mapOffDiagCoeffs
             {
                 FatalErrorIn
                 (
-                    "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
+                    "void extendedBlockLduMatrix<symmTensor>::"
+                    "mapOffDiagCoeffs\n"
+                    "(\n"
+                    "    const BlockLduMatrix<symmTensor>& blockLdum\n"
+                    ")"
                 )   << "Problem between ordinary block matrix and extended"
                     << " block matrix upper coeffs type morphing."
                     << abort(FatalError);
@@ -158,7 +167,11 @@ void Foam::extendedBlockLduMatrix<Foam::symmTensor>::mapOffDiagCoeffs
             {
                 FatalErrorIn
                 (
-                    "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
+                    "void extendedBlockLduMatrix<symmTensor>::"
+                    "mapOffDiagCoeffs\n"
+                    "(\n"
+                    "    const BlockLduMatrix<symmTensor>& blockLdum\n"
+                    ")"
                 )   << "Problem between ordinary block matrix and extended"
                     << " block matrix lower coeffs type morphing."
                     << abort(FatalError);
@@ -188,7 +201,8 @@ void Foam::extendedBlockLduMatrix<Foam::symmTensor>::mapOffDiagCoeffs
         if (upper.activeType() == blockCoeffBase::SCALAR)
         {
             // Helper type definition
-            typedef typename CoeffField<symmTensor>::scalarTypeField activeType;
+            typedef typename CoeffField<symmTensor>::scalarTypeField
+                activeType;
 
             // Get references to fields
             const activeType& activeUpper = upper.asScalar();
@@ -207,7 +221,8 @@ void Foam::extendedBlockLduMatrix<Foam::symmTensor>::mapOffDiagCoeffs
         else if (upper.activeType() == blockCoeffBase::LINEAR)
         {
             // Helper type definition
-            typedef typename CoeffField<symmTensor>::linearTypeField activeType;
+            typedef typename CoeffField<symmTensor>::linearTypeField
+                activeType;
 
             // Get references to fields
             const activeType& activeUpper = upper.asLinear();
@@ -227,7 +242,10 @@ void Foam::extendedBlockLduMatrix<Foam::symmTensor>::mapOffDiagCoeffs
         {
             FatalErrorIn
             (
-                "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
+                "void extendedBlockLduMatrix<symmTensor>::mapOffDiagCoeffs\n"
+                "(\n"
+                "    const BlockLduMatrix<symmTensor>& blockLdum\n"
+                ")"
             )   << "Problem between ordinary block matrix and extended"
                 << " block matrix upper/lower coeffs type morphing."
                 << abort(FatalError);
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/tensorExtendedBlockLduMatrix.C b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/tensorExtendedBlockLduMatrix.C
index 25efd56ea..dd7e768df 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/tensorExtendedBlockLduMatrix.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/extendedBlockLduMatrix/tensorExtendedBlockLduMatrix.C
@@ -41,8 +41,13 @@ void Foam::extendedBlockLduMatrix<Foam::tensor>::mapOffDiagCoeffs
 {
     if (blockLdum.diagonal())
     {
-        WarningIn("extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)")
-            << "Attempted to create extended lower/upper coeffs for block "
+        WarningIn
+        (
+            "void extendedBlockLduMatrix<tensor>::mapOffDiagCoeffs\n"
+            "(\n"
+            "    const BlockLduMatrix<tensor>& blockLdum\n"
+            ")"
+        )   << "Attempted to create extended lower/upper coeffs for block "
             << "matrix that is diagonal."
             << nl << endl;
     }
@@ -67,7 +72,8 @@ void Foam::extendedBlockLduMatrix<Foam::tensor>::mapOffDiagCoeffs
             if (upper.activeType() == blockCoeffBase::SCALAR)
             {
                 // Helper type definition
-                typedef typename CoeffField<tensor>::scalarTypeField activeType;
+                typedef typename CoeffField<tensor>::scalarTypeField
+                    activeType;
 
                 // Get references to fields
                 const activeType& activeUpper = upper.asScalar();
@@ -83,7 +89,8 @@ void Foam::extendedBlockLduMatrix<Foam::tensor>::mapOffDiagCoeffs
             else if (upper.activeType() == blockCoeffBase::LINEAR)
             {
                 // Helper type definition
-                typedef typename CoeffField<tensor>::linearTypeField activeType;
+                typedef typename CoeffField<tensor>::linearTypeField
+                    activeType;
 
                 // Get references to fields
                 const activeType& activeUpper = upper.asLinear();
@@ -100,7 +107,10 @@ void Foam::extendedBlockLduMatrix<Foam::tensor>::mapOffDiagCoeffs
             {
                 FatalErrorIn
                 (
-                    "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
+                    "void extendedBlockLduMatrix<tensor>::mapOffDiagCoeffs\n"
+                    "(\n"
+                    "    const BlockLduMatrix<tensor>& blockLdum\n"
+                    ")"
                 )   << "Problem between ordinary block matrix and extended"
                     << " block matrix upper coeffs type morphing."
                     << abort(FatalError);
@@ -121,7 +131,8 @@ void Foam::extendedBlockLduMatrix<Foam::tensor>::mapOffDiagCoeffs
             if (lower.activeType() == blockCoeffBase::SCALAR)
             {
                 // Helper type definition
-                typedef typename CoeffField<tensor>::scalarTypeField activeType;
+                typedef typename CoeffField<tensor>::scalarTypeField
+                    activeType;
 
                 // Get references to fields
                 const activeType& activeLower = lower.asScalar();
@@ -137,7 +148,8 @@ void Foam::extendedBlockLduMatrix<Foam::tensor>::mapOffDiagCoeffs
             else if (lower.activeType() == blockCoeffBase::LINEAR)
             {
                 // Helper type definition
-                typedef typename CoeffField<tensor>::linearTypeField activeType;
+                typedef typename CoeffField<tensor>::linearTypeField
+                    activeType;
 
                 // Get references to fields
                 const activeType& activeLower = lower.asLinear();
@@ -154,7 +166,10 @@ void Foam::extendedBlockLduMatrix<Foam::tensor>::mapOffDiagCoeffs
             {
                 FatalErrorIn
                 (
-                    "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
+                    "void extendedBlockLduMatrix<tensor>::mapOffDiagCoeffs\n"
+                    "(\n"
+                    "    const BlockLduMatrix<tensor>& blockLdum\n"
+                    ")"
                 )   << "Problem between ordinary block matrix and extended"
                     << " block matrix lower coeffs type morphing."
                     << abort(FatalError);
@@ -223,7 +238,10 @@ void Foam::extendedBlockLduMatrix<Foam::tensor>::mapOffDiagCoeffs
         {
             FatalErrorIn
             (
-                "extendedBlockLduMatrix(lduMatrix&, label, polyMesh&)"
+                "void extendedBlockLduMatrix<tensor>::mapOffDiagCoeffs\n"
+                "(\n"
+                "    const BlockLduMatrix<tensor>& blockLdum\n"
+                ")"
             )   << "Problem between ordinary block matrix and extended"
                 << " block matrix upper/lower coeffs type morphing."
                 << abort(FatalError);
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.C
index fff4a51c7..43f5a14be 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.C
@@ -67,7 +67,6 @@ void Foam::BlockGaussSeidelPrecon<Type>::calcInvDiag()
         // 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();
@@ -83,7 +82,7 @@ void Foam::BlockGaussSeidelPrecon<Type>::calcInvDiag()
         // Expand diagonal only to full square type and store into luDiag
         expandLinear(luDiag, lf);
 
-        // Keep only off-diagonal in ldDiag.
+        // Keep only off-diagonal in luDiag.
         // Note change of sign to avoid multiplication with -1 when moving
         // to the other side.  HJ, 20/Aug/2015
         luDiag -= activeDiag;
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.C
index 1f3589573..be722ae8d 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.C
@@ -65,8 +65,10 @@ void Foam::BlockILUCpPrecon<Type>::calcActiveTypeFactorization
         const extendedLduAddressing& addr = extBlockMatrix_.extendedLduAddr();
 
         // Get upper/lower extended addressing
-        const label* const __restrict__ uPtr = addr.extendedUpperAddr().begin();
-        const label* const __restrict__ lPtr = addr.extendedLowerAddr().begin();
+        const label* const __restrict__ uPtr =
+            addr.extendedUpperAddr().begin();
+        const label* const __restrict__ lPtr =
+            addr.extendedLowerAddr().begin();
 
         // Get extended owner start addressing
         const label* const __restrict__ ownStartPtr =
@@ -87,8 +89,8 @@ void Foam::BlockILUCpPrecon<Type>::calcActiveTypeFactorization
         LDUType* __restrict__ zPtr = z.begin();
         LDUType* __restrict__ wPtr = w.begin();
 
-        // Define start and end face ("virtual" face when extended addressing is
-        // used) of this row/column.
+        // Define start and end face ("virtual" face when extended addressing
+        // is used) of this row/column.
         register label fStart, fEnd, fLsrStart, fLsrEnd;
 
         // Crout LU factorization
@@ -505,17 +507,15 @@ Foam::BlockILUCpPrecon<Type>::BlockILUCpPrecon
 )
 :
     BlockLduPrecon<Type>(matrix),
-    preconDiag_(matrix.diag()),
-    p_(readLabel(dict.lookup("fillInLevel"))),
     extBlockMatrix_
     (
         matrix,
-        p_,
-        matrix.mesh().thisDb().time().template lookupObject<polyMesh>
+        matrix.mesh().lduAddr().extendedAddr
         (
-            polyMesh::defaultRegion
+            readLabel(dict.lookup("fillInLevel"))
         )
-    )
+    ),
+    preconDiag_(matrix.diag())
 {
     calcFactorization();
 }
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.H b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.H
index dc020ebbd..cf5c293b7 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.H
+++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.H
@@ -62,15 +62,12 @@ class BlockILUCpPrecon
 {
     // Private Data
 
+        //- Extended ldu matrix
+        mutable extendedBlockLduMatrix<Type> extBlockMatrix_;
+
         //- Preconditioned diagonal
         mutable CoeffField<Type> preconDiag_;
 
-        //- Fill in level
-        const label p_;
-
-        //- Extended lduMarix
-        mutable extendedBlockLduMatrix<Type> extBlockMatrix_;
-
 
     // Private Member Functions
 
diff --git a/src/foam/matrices/lduMatrix/lduAddressing/extendedLduAddressing/extendedLduAddressing.C b/src/foam/matrices/lduMatrix/lduAddressing/extendedLduAddressing/extendedLduAddressing.C
index 70f179a44..8e0aa6a37 100644
--- a/src/foam/matrices/lduMatrix/lduAddressing/extendedLduAddressing/extendedLduAddressing.C
+++ b/src/foam/matrices/lduMatrix/lduAddressing/extendedLduAddressing/extendedLduAddressing.C
@@ -25,6 +25,7 @@ License
 
 #include "extendedLduAddressing.H"
 #include "mapPolyMesh.H"
+#include "demandDrivenData.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -38,21 +39,10 @@ namespace Foam
 
 Foam::extendedLduAddressing::extendedLduAddressing
 (
-    const polyMesh& mesh,
     const lduAddressing& lduAddr,
     const label extensionLevel
 )
 :
-    MeshObject<polyMesh, extendedLduAddressing>(mesh),
-    lowerAddr_
-    (
-        labelList::subList
-        (
-            mesh.faceOwner(),
-            mesh.nInternalFaces()
-        )
-    ),
-    upperAddr_(mesh.faceNeighbour()),
     lduAddr_(lduAddr),
     p_(extensionLevel),
     extendedLowerPtr_(NULL),
@@ -146,6 +136,44 @@ void Foam::extendedLduAddressing::markNeighbours
 }
 
 
+void Foam::extendedLduAddressing::calcCellCells
+(
+    labelListList& cellCellAddr
+) const
+{
+    // Algorithm taken from primiteMeshCellCells.C in calcCellCells() member
+    // function.
+    labelList ncc(lduAddr_.size(), 0);
+
+    const unallocLabelList& own = lowerAddr();
+    const unallocLabelList& nei = upperAddr();
+
+    forAll (nei, faceI)
+    {
+        ncc[own[faceI]]++;
+        ncc[nei[faceI]]++;
+    }
+
+    // Size the list
+    cellCellAddr.setSize(ncc.size());
+
+    forAll (cellCellAddr, cellI)
+    {
+        cellCellAddr[cellI].setSize(ncc[cellI]);
+    }
+    ncc = 0;
+
+    forAll (nei, faceI)
+    {
+        label ownCellI = own[faceI];
+        label neiCellI = nei[faceI];
+
+        cellCellAddr[ownCellI][ncc[ownCellI]++] = neiCellI;
+        cellCellAddr[neiCellI][ncc[neiCellI]++] = ownCellI;
+    }
+}
+
+
 void Foam::extendedLduAddressing::calcExtendedLowerUpper() const
 {
     if (extendedLowerPtr_ && extendedUpperPtr_)
@@ -168,19 +196,16 @@ void Foam::extendedLduAddressing::calcExtendedLowerUpper() const
             << abort(FatalError);
     }
 
-    // Get polyMesh reference
-    const polyMesh& mesh = this->mesh();
-
-    // Allocate dynamic lists for extended owner/neighbour, which are later used
-    // to define ordinary labelLists (extendedLower, extendedUpper)
+    // Allocate dynamic lists for extended owner/neighbour, which are later
+    // used to define ordinary labelLists (extendedLower, extendedUpper)
     // Helper type definition
     typedef DynamicList<label> DynamicLabelList;
 
     // Get the number of faces
-    const label nFaces = upperAddr_.size();;
+    const label nFaces = upperAddr().size();
 
-    // Allocate extended owner/neighbour with 4*nFaces*extensionLevel as a guess
-    // (based on hex mesh assumption) to prevent excessive resizing
+    // Allocate extended owner/neighbour with 4*nFaces*extensionLevel as a
+    // guess (based on hex mesh assumption) to prevent excessive resizing
     DynamicLabelList dExtOwn(4*nFaces*p_);
     DynamicLabelList dExtNei(4*nFaces*p_);
 
@@ -188,10 +213,12 @@ void Foam::extendedLduAddressing::calcExtendedLowerUpper() const
     HashSet<label> extCellNbrs(32*p_);
 
     // Get number of cells in the mesh
-    const label nCells = mesh.nCells();
+    const label nCells = lduAddr_.size();
 
-    // Get a list of cells neighbouring each cell
-    const labelListList& cellCells = mesh.cellCells();
+    // Get a list of cells neighbouring each cell.
+    // Note optimisation to avoid a copy
+    labelListList cellCells;
+    calcCellCells(cellCells);
 
     // Loop through cells
     for (label cellI = 0; cellI < nCells; ++cellI)
@@ -234,8 +261,8 @@ void Foam::extendedLduAddressing::calcFaceMap() const
     }
 
     // Get reference to ordinary owner/neighbour addressing
-    const unallocLabelList& own = lduAddr().lowerAddr();
-    const unallocLabelList& nbr = lduAddr().upperAddr();
+    const unallocLabelList& own = lowerAddr();
+    const unallocLabelList& nbr = upperAddr();
 
     // Allocate memory for faceMap
     faceMapPtr_ = new labelList(own.size(), -1);
@@ -289,7 +316,7 @@ void Foam::extendedLduAddressing::calcExtendedLosort() const
     // Scan the extended neighbour list to find out how many times the cell
     // appears as a neighbour of the face. Done this way to avoid guessing
     // and resizing list
-    const label matrixSize = mesh().nCells();
+    const label matrixSize = lduAddr_.size();
     labelList nNbrOfFace(matrixSize, 0);
 
     const unallocLabelList& nbr = extendedUpperAddr();
@@ -429,7 +456,7 @@ void Foam::extendedLduAddressing::calcExtendedLosortStart() const
     }
 
     // Set up last lookup by hand
-    lsrtStart[this->mesh().nCells()] = nbr.size();
+    lsrtStart[lduAddr_.size()] = nbr.size();
 }
 
 
diff --git a/src/foam/matrices/lduMatrix/lduAddressing/extendedLduAddressing/extendedLduAddressing.H b/src/foam/matrices/lduMatrix/lduAddressing/extendedLduAddressing/extendedLduAddressing.H
index 83988b4ee..e6b1bf539 100644
--- a/src/foam/matrices/lduMatrix/lduAddressing/extendedLduAddressing/extendedLduAddressing.H
+++ b/src/foam/matrices/lduMatrix/lduAddressing/extendedLduAddressing/extendedLduAddressing.H
@@ -25,15 +25,15 @@ Class
     Foam::extendedLduAddressing
 
 Description
-    The class contains data and calculates demand driven extended ldu addressing
-    of arbitrary level. Level 0 is not allowed as it is the same as ordinary
-    lduAddressing. Level 1 considers neighbours of neighbours for a given cell,
-    etc.
+    The class contains data and calculates demand driven extended ldu
+    addressing of arbitrary level. Level 0 is not allowed as it is the same
+    as ordinary lduAddressing. Level 1 considers neighbours of neighbours
+    for a given cell, etc.
 
     The addressing is completely defined with extended variants of:
-    owner/neighbour, losort, ownerStart and losortStart: all being demand-driven
-    data. In order to define extended owner/neighbour lists, additional
-    "virtual" faces are assumed between cells not sharing a face.
+    owner/neighbour, losort, ownerStart and losortStart: all being
+    demand-driven data. In order to define extended owner/neighbour lists,
+    additional "virtual" faces are assumed between cells not sharing a face.
 
 Author
     Vuko Vukcevic, FMENA Zagreb. All rights reserved
@@ -46,9 +46,8 @@ SourceFiles
 #ifndef extendedLduAddressing_H
 #define extendedLduAddressing_H
 
-#include "MeshObject.H"
-#include "polyMesh.H"
 #include "lduAddressing.H"
+#include "HashSet.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -62,17 +61,9 @@ class mapPolyMesh;
 \*---------------------------------------------------------------------------*/
 
 class extendedLduAddressing
-:
-    public MeshObject<polyMesh, extendedLduAddressing>
 {
     // Private data
 
-        //- Lower as a subList of allOwner
-        labelList::subList lowerAddr_;
-
-        //- Upper as a reference to neighbour
-        const labelList& upperAddr_;
-
         //- Reference to basic lduAddressing
         const lduAddressing& lduAddr_;
 
@@ -121,7 +112,7 @@ class extendedLduAddressing
             //  Parameters:
             //  masterCellI - searching for extended neighbours of this cell
             //  nbrCells - neighbouring cells of a cell (based on level)
-            //  cellCells - mesh data to avoid function calls
+            //  cellCells - passing as a parameter for minor optimisation
             //  extendedNbrs - HashSet containing labels of extended neighbours
             //  curLevel - current level needed for recursion control
             void markNeighbours
@@ -136,6 +127,10 @@ class extendedLduAddressing
 
         // Calculation of demand driven data
 
+            //- Calculate cells cells using lower/upper addressing
+            //  into provided storage
+            void calcCellCells(labelListList& cc) const;
+
             //- Calculate extended lower/upper
             void calcExtendedLowerUpper() const;
 
@@ -160,10 +155,9 @@ public:
 
     // Constructors
 
-        //- Construct given polyMesh
+        //- Construct given lduAddressing
         extendedLduAddressing
         (
-            const polyMesh& mesh,
             const lduAddressing& lduAddr,
             const label extensionLevel
         );
@@ -178,22 +172,22 @@ public:
 
         // Access
 
+            //- Return basic lduAddressing
+            const lduAddressing& lduAddr() const
+            {
+                return lduAddr_;
+            }
+
             //- Return lower addressing (i.e. lower label = upper triangle)
             const unallocLabelList& lowerAddr() const
             {
-                return lowerAddr_;
+                return lduAddr_.lowerAddr();
             }
 
             //- Return upper addressing (i.e. upper label)
             const unallocLabelList& upperAddr() const
             {
-                return upperAddr_;
-            }
-
-            //- Return basic lduAddressing
-            const lduAddressing& lduAddr() const
-            {
-                return lduAddr_;
+                return lduAddr_.upperAddr();
             }
 
             //- Return extended lower addressing
diff --git a/src/foam/matrices/lduMatrix/lduAddressing/lduAddressing.C b/src/foam/matrices/lduMatrix/lduAddressing/lduAddressing.C
index 28f627ce8..5e975409b 100644
--- a/src/foam/matrices/lduMatrix/lduAddressing/lduAddressing.C
+++ b/src/foam/matrices/lduMatrix/lduAddressing/lduAddressing.C
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "lduAddressing.H"
+#include "extendedLduAddressing.H"
 #include "demandDrivenData.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
@@ -169,6 +170,18 @@ void Foam::lduAddressing::calcLosortStart() const
 }
 
 
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::lduAddressing::lduAddressing(const label nEqns)
+:
+    size_(nEqns),
+    losortPtr_(NULL),
+    ownerStartPtr_(NULL),
+    losortStartPtr_(NULL),
+    extendedAddr_(5)
+{}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::lduAddressing::~lduAddressing()
@@ -248,4 +261,27 @@ Foam::label Foam::lduAddressing::triIndex(const label a, const label b) const
 }
 
 
+const Foam::extendedLduAddressing&
+Foam::lduAddressing::extendedAddr(const label p) const
+{
+    if (p == 0 || p > 4)
+    {
+        FatalErrorIn
+        (
+            "const Foam::extendedLduAddressing& "
+            "Foam::lduAddressing::extendedAddr(const label p) const"
+        )   << "Currently supported extended addressing fill-in only "
+            << "between order 1 and 4"
+            << abort(FatalError);
+    }
+
+    if (!extendedAddr_.set(p))
+    {
+        extendedAddr_.set(p, new extendedLduAddressing(*this, p));
+    }
+
+    return extendedAddr_[p];
+}
+
+
 // ************************************************************************* //
diff --git a/src/foam/matrices/lduMatrix/lduAddressing/lduAddressing.H b/src/foam/matrices/lduMatrix/lduAddressing/lduAddressing.H
index f0818c478..709fa9caf 100644
--- a/src/foam/matrices/lduMatrix/lduAddressing/lduAddressing.H
+++ b/src/foam/matrices/lduMatrix/lduAddressing/lduAddressing.H
@@ -104,8 +104,11 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
+class extendedLduAddressing;
+
 /*---------------------------------------------------------------------------*\
-                           Class lduAddressing Declaration
+                       Class lduAddressing Declaration
 \*---------------------------------------------------------------------------*/
 
 class lduAddressing
@@ -127,6 +130,9 @@ class lduAddressing
         //- Losort start addressing
         mutable labelList* losortStartPtr_;
 
+        //- Extended addressing with p-order fill-in
+        mutable PtrList<extendedLduAddressing> extendedAddr_;
+
 
     // Private Member Functions
 
@@ -149,18 +155,11 @@ class lduAddressing
 public:
 
     // Constructor
-    lduAddressing(const label nEqns)
-    :
-        size_(nEqns),
-        losortPtr_(NULL),
-        ownerStartPtr_(NULL),
-        losortStartPtr_(NULL)
-    {}
+    lduAddressing(const label nEqns);
 
 
-    // Destructor
-
-        virtual ~lduAddressing();
+    //- Destructor
+    virtual ~lduAddressing();
 
 
     // Member Functions
@@ -200,6 +199,9 @@ public:
 
         //- Return off-diagonal index given owner and neighbour label
         label triIndex(const label a, const label b) const;
+
+        //- Return extended addressing given p index
+        const extendedLduAddressing& extendedAddr(const label p) const;
 };
 
 
diff --git a/src/foam/matrices/lduMatrix/lduMatrix/extendedLduMatrix/extendedLduMatrix.C b/src/foam/matrices/lduMatrix/lduMatrix/extendedLduMatrix/extendedLduMatrix.C
index a7be169e5..711a2f28d 100644
--- a/src/foam/matrices/lduMatrix/lduMatrix/extendedLduMatrix/extendedLduMatrix.C
+++ b/src/foam/matrices/lduMatrix/lduMatrix/extendedLduMatrix/extendedLduMatrix.C
@@ -38,20 +38,11 @@ namespace Foam
 Foam::extendedLduMatrix::extendedLduMatrix
 (
     const lduMatrix& ldum,
-    const label extensionLevel,
-    const polyMesh& polyMesh
+    const extendedLduAddressing& extLduAddr
 )
 :
     basicLduMatrix_(ldum),
-    extLduAddr_
-    (
-        extendedLduAddressing::New
-        (
-            polyMesh,
-            ldum.lduAddr(),
-            extensionLevel
-        )
-    ),
+    extLduAddr_(extLduAddr),
     extendedLowerPtr_(NULL),
     extendedUpperPtr_(NULL)
 {
@@ -102,10 +93,10 @@ Foam::extendedLduMatrix::extendedLduMatrix
         const label nExtFaces = extLduAddr_.extendedUpperAddr().size();
 
         // Allocate extended upper and lower
-        extendedUpperPtr_ = new scalarField(nExtFaces, 0.0);
+        extendedUpperPtr_ = new scalarField(nExtFaces, scalar(0));
         scalarField& extUpper = *extendedUpperPtr_;
 
-        extendedLowerPtr_ = new scalarField(nExtFaces, 0.0);
+        extendedLowerPtr_ = new scalarField(nExtFaces, scalar(0));
         scalarField& extLower = *extendedLowerPtr_;
 
         // Get upper and lower coeffs from underlying lduMatrix
@@ -154,7 +145,7 @@ Foam::scalarField& Foam::extendedLduMatrix::extendedLower()
             extendedLowerPtr_ = new scalarField
             (
                 extLduAddr_.extendedLowerAddr().size(),
-                0.0
+                scalar(0)
             );
         }
     }
@@ -176,7 +167,7 @@ Foam::scalarField& Foam::extendedLduMatrix::extendedUpper()
             extendedUpperPtr_ = new scalarField
             (
                 extLduAddr_.extendedUpperAddr().size(),
-                0.0
+                scalar(0)
             );
         }
     }
diff --git a/src/foam/matrices/lduMatrix/lduMatrix/extendedLduMatrix/extendedLduMatrix.H b/src/foam/matrices/lduMatrix/lduMatrix/extendedLduMatrix/extendedLduMatrix.H
index 7357f603e..d1300200a 100644
--- a/src/foam/matrices/lduMatrix/lduMatrix/extendedLduMatrix/extendedLduMatrix.H
+++ b/src/foam/matrices/lduMatrix/lduMatrix/extendedLduMatrix/extendedLduMatrix.H
@@ -39,9 +39,8 @@ SourceFiles
 #define extendedLduMatrix_H
 
 #include "extendedLduAddressing.H"
-#include "polyMesh.H"
-#include "className.H"
 #include "lduMatrix.H"
+#include "className.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -88,12 +87,15 @@ public:
     // Constructors
 
         //- Construct given lduMatrix, extension level and polyMesh
-        extendedLduMatrix(const lduMatrix&, const label, const polyMesh&);
+        extendedLduMatrix
+        (
+            const lduMatrix&,
+            const extendedLduAddressing&
+        );
 
 
-    // Destructor
-
-        ~extendedLduMatrix();
+    //- Destructor
+    ~extendedLduMatrix();
 
 
     // Member functions
diff --git a/src/foam/primitives/VectorN/TensorNI.H b/src/foam/primitives/VectorN/TensorNI.H
index 764024d1c..9e80eec06 100644
--- a/src/foam/primitives/VectorN/TensorNI.H
+++ b/src/foam/primitives/VectorN/TensorNI.H
@@ -104,13 +104,14 @@ inline TensorN<Cmpt, length> TensorN<Cmpt, length>::T() const
     return transpose;
 }
 
+
 //- Return tensor diagonal
 template <class Cmpt, int length>
 inline DiagTensorN<Cmpt, length> TensorN<Cmpt, length>::diag() const
 {
     DiagTensorN<Cmpt, length> dt;
 
-    label diagI=0;
+    label diagI = 0;
     for (label i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
     {
         dt[i] = this->v_[diagI];
@@ -120,6 +121,7 @@ inline DiagTensorN<Cmpt, length> TensorN<Cmpt, length>::diag() const
     return dt;
 }
 
+
 //- Negative sum the vertical off-diagonal components
 template <class Cmpt, int length>
 inline TensorN<Cmpt, length> TensorN<Cmpt, length>::negSumDiag() const
@@ -153,6 +155,7 @@ inline TensorN<Cmpt, length> TensorN<Cmpt, length>::negSumDiag() const
     return negsumdiag;
 }
 
+
 //- Assign to a SphericalTensorN
 template <class Cmpt, int length>
 inline void
@@ -502,7 +505,11 @@ operator+(const TensorN<Cmpt, length>& t1, const TensorN<Cmpt, length>& t2)
 //- Addition of TensorN and DiagTensorN
 template <class Cmpt, int length>
 inline TensorN<Cmpt, length>
-operator+(const TensorN<Cmpt, length>& t1, const DiagTensorN<Cmpt, length>& dt2)
+operator+
+(
+    const TensorN<Cmpt, length>& t1,
+    const DiagTensorN<Cmpt, length>& dt2
+)
 {
     TensorN<Cmpt, length> result(t1);
 
@@ -520,7 +527,11 @@ operator+(const TensorN<Cmpt, length>& t1, const DiagTensorN<Cmpt, length>& dt2)
 //- Addition of DiagTensorN and TensorN
 template <class Cmpt, int length>
 inline TensorN<Cmpt, length>
-operator+(const DiagTensorN<Cmpt, length>& dt1, const TensorN<Cmpt, length>& t2)
+operator+
+(
+    const DiagTensorN<Cmpt, length>& dt1,
+    const TensorN<Cmpt, length>& t2
+)
 {
     TensorN<Cmpt, length> result(t2);
 
@@ -602,7 +613,11 @@ operator-(const TensorN<Cmpt, length>& t1, const TensorN<Cmpt, length>& t2)
 //- Subtraction of TensorN and DiagTensorN
 template <class Cmpt, int length>
 inline TensorN<Cmpt, length>
-operator-(const TensorN<Cmpt, length>& t1, const DiagTensorN<Cmpt, length>& dt2)
+operator-
+(
+    const TensorN<Cmpt, length>& t1,
+    const DiagTensorN<Cmpt, length>& dt2
+)
 {
     TensorN<Cmpt, length> result(t1);
 
@@ -620,7 +635,11 @@ operator-(const TensorN<Cmpt, length>& t1, const DiagTensorN<Cmpt, length>& dt2)
 //- Subtraction of DiagTensorN and TensorN
 template <class Cmpt, int length>
 inline TensorN<Cmpt, length>
-operator-(const DiagTensorN<Cmpt, length>& dt1, const TensorN<Cmpt, length>& t2)
+operator-
+(
+    const DiagTensorN<Cmpt, length>& dt1,
+    const TensorN<Cmpt, length>& t2
+)
 {
     TensorN<Cmpt, length> result(-t2);
 
@@ -709,7 +728,11 @@ operator/(const TensorN<Cmpt, length>& t1, const TensorN<Cmpt, length>& t2)
 //- Inner Product of a DiagTensorN and an inverse TensorN
 template <class Cmpt, int length>
 inline TensorN<Cmpt, length>
-operator/(const DiagTensorN<Cmpt, length>& dt1, const TensorN<Cmpt, length>& t2)
+operator/
+(
+    const DiagTensorN<Cmpt, length>& dt1,
+    const TensorN<Cmpt, length>& t2
+)
 {
     return dt1 & inv(t2);
 }
@@ -718,7 +741,11 @@ operator/(const DiagTensorN<Cmpt, length>& dt1, const TensorN<Cmpt, length>& t2)
 //- Inner Product of a TensorN and an inverse DiagTensorN
 template <class Cmpt, int length>
 inline TensorN<Cmpt, length>
-operator/(const TensorN<Cmpt, length>& t1, const DiagTensorN<Cmpt, length>& dt2)
+operator/
+(
+    const TensorN<Cmpt, length>& t1,
+    const DiagTensorN<Cmpt, length>& dt2
+)
 {
     return t1 & inv(dt2);
 }
@@ -972,6 +999,7 @@ inline DiagTensorN<Cmpt, length> diag(const TensorN<Cmpt, length>& t)
     return t.diag();
 }
 
+
 //- Return tensor diagonal
 template <class Cmpt, int length>
 inline TensorN<Cmpt, length> negSumDiag(const TensorN<Cmpt, length>& t)
diff --git a/src/foam/primitives/VectorN/VectorN.H b/src/foam/primitives/VectorN/VectorN.H
index a4375a90b..dfca611cc 100644
--- a/src/foam/primitives/VectorN/VectorN.H
+++ b/src/foam/primitives/VectorN/VectorN.H
@@ -70,6 +70,8 @@ public:
         static const char* const typeName;
         static const VectorN zero;
         static const VectorN one;
+        static const VectorN max;
+        static const VectorN min;
 
 
     // Constructors
diff --git a/src/foam/primitives/VectorN/VectorNI.H b/src/foam/primitives/VectorN/VectorNI.H
index 2dfd8fcf4..a0759b0eb 100644
--- a/src/foam/primitives/VectorN/VectorNI.H
+++ b/src/foam/primitives/VectorN/VectorNI.H
@@ -43,6 +43,12 @@ const VectorN<Cmpt, length> VectorN<Cmpt, length>::zero(0);
 template <class Cmpt, int length>
 const VectorN<Cmpt, length> VectorN<Cmpt, length>::one(1);
 
+template <class Cmpt, int length>
+const VectorN<Cmpt, length> VectorN<Cmpt, length>::max(VGREAT);
+
+template <class Cmpt, int length>
+const VectorN<Cmpt, length> VectorN<Cmpt, length>::min(-VGREAT);
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
diff --git a/src/foam/primitives/VectorN/expandContract/ExpandTensorN.H b/src/foam/primitives/VectorN/expandContract/ExpandTensorN.H
index d78304882..60364aab9 100644
--- a/src/foam/primitives/VectorN/expandContract/ExpandTensorN.H
+++ b/src/foam/primitives/VectorN/expandContract/ExpandTensorN.H
@@ -75,6 +75,9 @@ UNARY_TEMPLATE_FUNCTION(TensorN, VectorN, expandLinear)
 UNARY_TEMPLATE_FUNCTION(DiagTensorN, VectorN, expandLinear)
 UNARY_TEMPLATE_FUNCTION(SphericalTensorN, VectorN, expandLinear)
 
+UNARY_TEMPLATE_FUNCTION(VectorN, TensorN, sumToDiag)
+UNARY_TEMPLATE_FUNCTION(VectorN, TensorN, sumMagToDiag)
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/foam/primitives/VectorN/expandContract/ExpandTensorNI.H b/src/foam/primitives/VectorN/expandContract/ExpandTensorNI.H
index 34dabb6e2..f36cc8e66 100644
--- a/src/foam/primitives/VectorN/expandContract/ExpandTensorNI.H
+++ b/src/foam/primitives/VectorN/expandContract/ExpandTensorNI.H
@@ -37,7 +37,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-//- Return the Average of a vector of as a scalar
+//- Return the average of a vector of as a scalar
 template <class Cmpt, int length>
 inline void contractScalar(Cmpt& result, const VectorN<Cmpt, length>& t)
 {
@@ -48,7 +48,7 @@ inline void contractScalar(Cmpt& result, const VectorN<Cmpt, length>& t)
         result += t[i];
     }
 
-//- Modified 2009/11/3 by I. Clifford
+    // Modified 2009/11/3 by I. Clifford
     result *= 1.0/VectorN<Cmpt, length>::nComponents;
 }
 
@@ -76,7 +76,7 @@ inline void contractScalar(Cmpt& result, const TensorN<Cmpt, length>& t)
         j += TensorN<Cmpt, length>::rowLength+1;
     }
 
-//- Modified 2009/11/3 by I. Clifford
+    // Modified 2009/11/3 by I. Clifford
     result *= 1.0/TensorN<Cmpt, length>::rowLength;
 }
 
@@ -118,7 +118,11 @@ inline Cmpt contractScalar(const DiagTensorN<Cmpt, length>& t)
 
 //- Return the diagonal of a SphericalTensorN as a scalar
 template <class Cmpt, int length>
-inline void contractScalar(Cmpt& result, const SphericalTensorN<Cmpt, length>& t)
+inline void contractScalar
+(
+    Cmpt& result,
+    const SphericalTensorN<Cmpt, length>& t
+)
 {
     result = t[0];
 }
@@ -134,7 +138,7 @@ inline Cmpt contractScalar(const SphericalTensorN<Cmpt, length>& t)
 }
 
 
-//- Return the diagonal of a TensorN as a vector
+//- Return the diagonal of a TensorN as a vectorN
 template <class Cmpt, int length>
 inline void contractLinear
 (
@@ -142,16 +146,16 @@ inline void contractLinear
     const TensorN<Cmpt, length>& t
 )
 {
-    int j=0;
+    int j = 0;
     for (int i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
     {
         result[i] = t[j];
-        j += TensorN<Cmpt, length>::rowLength+1;
+        j += TensorN<Cmpt, length>::rowLength + 1;
     }
 }
 
 
-//- Return the diagonal of a TensorN as a vector
+//- Return the diagonal of a TensorN as a vectorN
 template <class Cmpt, int length>
 inline VectorN<Cmpt, length> contractLinear(const TensorN<Cmpt, length>& t)
 {
@@ -161,7 +165,7 @@ inline VectorN<Cmpt, length> contractLinear(const TensorN<Cmpt, length>& t)
 }
 
 
-//- Return the diagonal of a DiagTensorN as a vector
+//- Return the diagonal of a DiagTensorN as a vectorN
 template <class Cmpt, int length>
 inline void contractLinear
 (
@@ -176,7 +180,7 @@ inline void contractLinear
 }
 
 
-//- Return the diagonal of a DiagTensorN as a vector
+//- Return the diagonal of a DiagTensorN as a vectorN
 template <class Cmpt, int length>
 inline VectorN<Cmpt, length> contractLinear(const DiagTensorN<Cmpt, length>& t)
 {
@@ -186,7 +190,7 @@ inline VectorN<Cmpt, length> contractLinear(const DiagTensorN<Cmpt, length>& t)
 }
 
 
-//- Return the diagonal of a SphericalTensorN as a vector
+//- Return the diagonal of a SphericalTensorN as a vectorN
 template <class Cmpt, int length>
 inline void contractLinear
 (
@@ -201,7 +205,7 @@ inline void contractLinear
 }
 
 
-//- Return the diagonal of a SphericalTensorN as a vector
+//- Return the diagonal of a SphericalTensorN as a vectorN
 template <class Cmpt, int length>
 inline VectorN<Cmpt, length> contractLinear
 (
@@ -311,6 +315,64 @@ inline void expandLinear
 }
 
 
+//- Sum row elements of a TensorN as a vectorN
+template <class Cmpt, int length>
+inline void sumToDiag
+(
+    VectorN<Cmpt, length>& result,
+    const TensorN<Cmpt, length>& t
+)
+{
+    result = VectorN<Cmpt, length>::zero;
+
+    for (int i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
+    {
+        for (int j = 0; j < TensorN<Cmpt, length>::rowLength; j++)
+        {
+            result[i] += t(i, j);
+        }
+    }
+}
+
+
+//- Sum row elements of a TensorN as a vectorN
+template <class Cmpt, int length>
+inline VectorN<Cmpt, length> sumToDiag(const TensorN<Cmpt, length>& t)
+{
+    VectorN<Cmpt, length> result;
+    sumToDiag(result, t);
+    return result;
+}
+
+
+//- Sum magnitudes of row elements of a TensorN as a vectorN
+template <class Cmpt, int length>
+inline void sumMagToDiag
+(
+    VectorN<Cmpt, length>& result,
+    const TensorN<Cmpt, length>& t
+)
+{
+    for (int i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
+    {
+        for (int j = 0; j < TensorN<Cmpt, length>::rowLength; j++)
+        {
+            result[i] += Foam::mag(t(i, j));
+        }
+    }
+}
+
+
+//- Sum magnitudes of row elements of a TensorN as a vectorN
+template <class Cmpt, int length>
+inline VectorN<Cmpt, length> sumMagToDiag(const TensorN<Cmpt, length>& t)
+{
+    VectorN<Cmpt, length> result;
+    sumMagToDiag(result, t);
+    return result;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/foam/primitives/expandContract/expandTensor.H b/src/foam/primitives/expandContract/expandTensor.H
index 525643f88..ff68c647f 100644
--- a/src/foam/primitives/expandContract/expandTensor.H
+++ b/src/foam/primitives/expandContract/expandTensor.H
@@ -110,6 +110,54 @@ inline void expandLinear(Tensor<Cmpt>& result, const Vector<Cmpt>& v)
 }
 
 
+//- Sum row elements of a tensor as a vector
+template <class Cmpt>
+inline void sumToDiag(Vector<Cmpt>& result, const Tensor<Cmpt>& t)
+{
+    result = Vector<Cmpt>
+    (
+        t.xx() + t.xy() + t.xz(),
+        t.yx() + t.yy() + t.yz(),
+        t.zx() + t.zy() + t.zz()
+    );
+}
+
+
+//- Sum row elements of a TensorN as a vectorN
+template <class Cmpt>
+inline Vector<Cmpt> sumToDiag(const Tensor<Cmpt>& t)
+{
+    Vector<Cmpt> result;
+    sumToDiag(result, t);
+    return result;
+}
+
+
+//- Sum row elements of a tensor as a vector
+template <class Cmpt>
+inline void sumMagToDiag(Vector<Cmpt>& result, const Tensor<Cmpt>& t)
+{
+    result = Vector<Cmpt>
+    (
+        Foam::mag(t.xx()) + Foam::mag(t.xy()) + Foam::mag(t.xz()),
+        Foam::mag(t.yx()) + Foam::mag(t.yy()) + Foam::mag(t.yz()),
+        Foam::mag(t.zx()) + Foam::mag(t.zy()) + Foam::mag(t.zz())
+    );
+}
+
+
+//- Sum magnitudes of row elements of a TensorN as a vectorN
+template <class Cmpt>
+inline Vector<Cmpt> sumMagToDiag(const Tensor<Cmpt>& t)
+{
+    Vector<Cmpt> result;
+    sumMagToDiag(result, t);
+    return result;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/lduSolvers/lduPrecon/ILUCp/ILUCp.C b/src/lduSolvers/lduPrecon/ILUCp/ILUCp.C
index 969cc0896..47eebdf4f 100644
--- a/src/lduSolvers/lduPrecon/ILUCp/ILUCp.C
+++ b/src/lduSolvers/lduPrecon/ILUCp/ILUCp.C
@@ -50,8 +50,8 @@ namespace Foam
         addILUCpPreconditionerAsymMatrixConstructorToTable_;
 
     // Add to symmetric constructor table as well. Cholesky with fill in would
-    // yield the same sparseness pattern as the original matrix, hence it is not
-    // implemented. VV, 10/Sep/2015.
+    // yield the same sparseness pattern as the original matrix,
+    // hence it is not implemented. VV, 10/Sep/2015.
     lduPreconditioner::
         addsymMatrixConstructorToTable<ILUCp>
         addILUCpPreconditionerSymMatrixConstructorToTable_;
@@ -69,8 +69,10 @@ void Foam::ILUCp::calcFactorization()
         const extendedLduAddressing& addr = extMatrix_.extendedLduAddr();
 
         // Get upper/lower extended addressing
-        const label* const __restrict__ uPtr = addr.extendedUpperAddr().begin();
-        const label* const __restrict__ lPtr = addr.extendedLowerAddr().begin();
+        const label* const __restrict__ uPtr =
+            addr.extendedUpperAddr().begin();
+        const label* const __restrict__ lPtr =
+            addr.extendedLowerAddr().begin();
 
         // Get extended owner start addressing
         const label* const __restrict__ ownStartPtr =
@@ -94,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
-        // used) of this row/column.
+        // Define start and end face ("virtual" face when extended addressing
+        // is used) of this row/column.
         register label fStart, fEnd, fLsrStart, fLsrEnd;
 
         // Crout LU factorization
@@ -239,17 +241,15 @@ Foam::ILUCp::ILUCp
         coupleIntCoeffs,
         interfaces
     ),
-    preconDiag_(matrix_.diag()),
-    p_(readLabel(dict.lookup("fillInLevel"))),
     extMatrix_
     (
         matrix,
-        p_,
-        matrix.mesh().thisDb().time().lookupObject
-        <
-            polyMesh
-        >(polyMesh::defaultRegion)
+        matrix.mesh().lduAddr().extendedAddr
+        (
+            readLabel(dict.lookup("fillInLevel"))
+        )
     ),
+    preconDiag_(matrix_.diag()),
     zDiag_(0),
     z_(preconDiag_.size(), 0),
     w_(preconDiag_.size(), 0)
diff --git a/src/lduSolvers/lduPrecon/ILUCp/ILUCp.H b/src/lduSolvers/lduPrecon/ILUCp/ILUCp.H
index 274118b3c..6e671d356 100644
--- a/src/lduSolvers/lduPrecon/ILUCp/ILUCp.H
+++ b/src/lduSolvers/lduPrecon/ILUCp/ILUCp.H
@@ -51,7 +51,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                       Class ILUCp Declaration
+                           Class ILUCp Declaration
 \*---------------------------------------------------------------------------*/
 
 class ILUCp
@@ -60,15 +60,12 @@ class ILUCp
 {
     // Private Data
 
+        //- Extended ldu matrix
+        extendedLduMatrix extMatrix_;
+
         //- Preconditioned diagonal
         scalarField preconDiag_;
 
-        //- Fill in level
-        const label p_;
-
-        //- Extended lduMatrix
-        extendedLduMatrix extMatrix_;
-
         //- Temporary working diagonal
         scalar zDiag_;