Merge remote branch 'remotes/origin/blockCoupledFVM'

src/OpenFOAM/Make/files

@@ -578,4 +578,55 @@ $(writers)/gnuplotGraph/gnuplotGraph.C
 $(writers)/xmgrGraph/xmgrGraph.C
 $(writers)/jplotGraph/jplotGraph.C
 
+//- Block Matrix Additions
+//- I. Clifford 03-12-2009
+primitives/BlockCoeff/blockCoeffBase.C
+primitives/BlockCoeff/scalarBlockCoeff.C
+primitives/BlockCoeff/sphericalTensorBlockCoeff.C
+primitives/BlockCoeff/symmTensorBlockCoeff.C
+primitives/BlockCoeff/tensorBlockCoeff.C
+
+fields/expandContract/expandTensorField.C
+
+fields/CoeffField/scalarCoeffField.C
+// fields/CoeffField/sphericalTensorCoeffField.C
+// fields/CoeffField/symmTensorCoeffField.C
+fields/CoeffField/tensorCoeffField.C
+
+matrices/blockLduMatrix/BlockLduMatrix/blockLduMatrices.C
+matrices/blockLduMatrix/BlockLduMatrix/scalarBlockLduMatrix.C
+// matrices/blockLduMatrix/BlockLduMatrix/sphericalTensorBlockLduMatrix.C
+// matrices/blockLduMatrix/BlockLduMatrix/symmTensorBlockLduMatrix.C
+matrices/blockLduMatrix/BlockLduMatrix/tensorBlockLduMatrix.C
+matrices/blockLduMatrix/BlockLduMatrix/BlockConstraint/scalarBlockConstraint.C
+
+matrices/blockLduMatrix/BlockLduPrecons/BlockLduPrecon/blockLduPrecons.C
+matrices/blockLduMatrix/BlockLduPrecons/BlockNoPrecon/blockNoPrecons.C
+
+matrices/blockLduMatrix/BlockLduPrecons/BlockDiagonalPrecon/scalarBlockDiagonalPrecon.C
+matrices/blockLduMatrix/BlockLduPrecons/BlockDiagonalPrecon/tensorBlockDiagonalPrecon.C
+matrices/blockLduMatrix/BlockLduPrecons/BlockDiagonalPrecon/blockDiagonalPrecons.C
+
+matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/scalarBlockGaussSeidelPrecon.C
+matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/tensorBlockGaussSeidelPrecon.C
+matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/blockGaussSeidelPrecons.C
+
+matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/scalarBlockCholeskyPrecon.C
+matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/tensorBlockCholeskyPrecon.C
+matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/blockCholeskyPrecons.C
+
+matrices/blockLduMatrix/BlockLduSmoothers/BlockLduSmoother/blockLduSmoothers.C
+// matrices/blockLduMatrix/BlockLduSmoothers/BlockGaussSeidelSmoother/blockGaussSeidelSmoothers.C
+matrices/blockLduMatrix/BlockLduSmoothers/BlockILUSmoother/blockILUSmoothers.C
+
+matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/blockLduSolvers.C
+matrices/blockLduMatrix/BlockLduSolvers/BlockDiagonal/blockDiagonalSolvers.C
+// matrices/blockLduMatrix/BlockLduSolvers/Segregated/segregatedSolvers.C
+matrices/blockLduMatrix/BlockLduSolvers/BlockGaussSeidel/blockGaussSeidelSolvers.C
+matrices/blockLduMatrix/BlockLduSolvers/BlockCG/blockCGSolvers.C
+matrices/blockLduMatrix/BlockLduSolvers/BlockBiCGStab/blockBiCGStabSolvers.C
+matrices/blockLduMatrix/BlockLduSolvers/BlockGMRES/blockGMRESSolvers.C
+
+matrices/blockLduMatrix/BlockLduMatrix/BlockLduInterface/BlockLduInterfaceFields.C
+
 LIB = $(FOAM_LIBBIN)/libOpenFOAM

src/OpenFOAM/fields/CoeffField/CoeffField.C

@@ -998,20 +998,234 @@ void Foam::CoeffField<Type>::addSubset
         }
     }
     else
+    {
+        if (f.activeType() == blockCoeffBase::SQUARE)
+        {
+            squareTypeField& localF = this->toSquare();
+
+            const squareTypeField& ff = f.asSquare();
+
+            forAll (ff, ffI)
+            {
+                localF[addr[ffI]] += ff[ffI];
+            }
+        }
+        if (f.activeType() == blockCoeffBase::LINEAR)
+        {
+            linearTypeField& localF = this->toLinear();
+
+            const linearTypeField& ff = f.asLinear();
+
+            forAll (ff, ffI)
+            {
+                localF[addr[ffI]] += ff[ffI];
+            }
+        }
+        else if (f.activeType() == blockCoeffBase::SCALAR)
+        {
+            scalarTypeField& localF = this->toScalar();
+
+            const scalarTypeField& ff = f.asScalar();
+
+            forAll (ff, ffI)
+            {
+                localF[addr[ffI]] += ff[ffI];
+            }
+        }
+//         FatalErrorIn
+//         (
+//             "template<class Type>\n"
+//             "void Foam::CoeffField<Type>::addSubset\n"
+//             "(\n"
+//             "    const CoeffField<Type>& f,\n"
+//             "    const labelList addr\n"
+//             ")"
+//         )   << "Incompatible combination of types"
+//             << abort(FatalError);
+    }
+}
+
+
+template<class Type>
+void Foam::CoeffField<Type>::subtractSubset
+(
+    const CoeffField<Type>& f,
+    const labelList& addr
+)
+{
+    // Check sizes
+    if (f.size() != addr.size())
     {
         FatalErrorIn
         (
             "template<class Type>\n"
-            "void Foam::CoeffField<Type>::addSubset\n"
+            "void Foam::CoeffField<Type>::subtractSubset\n"
             "(\n"
             "    const CoeffField<Type>& f,\n"
             "    const labelList addr\n"
             ")"
-        )   << "Incompatible combination of types"
+        )   << "Incompatible sizes: " << f.size() << " and " << addr.size()
             << abort(FatalError);
     }
-}
 
+    if (this->activeType() == blockCoeffBase::SQUARE)
+    {
+        if (f.activeType() == blockCoeffBase::SQUARE)
+        {
+            squareTypeField& localF = this->asSquare();
+
+            const squareTypeField& ff = f.asSquare();
+
+            forAll (ff, ffI)
+            {
+                localF[addr[ffI]] -= ff[ffI];
+            }
+        }
+        else if (f.activeType() == blockCoeffBase::LINEAR)
+        {
+            squareTypeField& localF = this->asSquare();
+
+            squareTypeField ff(f.size());
+            expandLinear(ff, f.asLinear());
+
+            forAll (ff, ffI)
+            {
+                localF[addr[ffI]] -= ff[ffI];
+            }
+        }
+        else if (f.activeType() == blockCoeffBase::SCALAR)
+        {
+            squareTypeField& localF = this->asSquare();
+
+            squareTypeField ff(f.size());
+            expandScalar(ff, f.asScalar());
+
+            forAll (ff, ffI)
+            {
+                localF[addr[ffI]] -= ff[ffI];
+            }
+        }
+    }
+    else if (this->activeType() == blockCoeffBase::LINEAR)
+    {
+        if (f.activeType() == blockCoeffBase::SQUARE)
+        {
+            squareTypeField& localF = this->asSquare();
+
+            const squareTypeField& ff = f.asSquare();
+
+            forAll (ff, ffI)
+            {
+                localF[addr[ffI]] -= ff[ffI];
+            }
+        }
+        if (f.activeType() == blockCoeffBase::LINEAR)
+        {
+            linearTypeField& localF = this->asLinear();
+
+            const linearTypeField& ff = f.asLinear();
+
+            forAll (ff, ffI)
+            {
+                localF[addr[ffI]] -= ff[ffI];
+            }
+        }
+        else if (f.activeType() == blockCoeffBase::SCALAR)
+        {
+            linearTypeField& localF = this->asLinear();
+
+            linearTypeField ff(f.size());
+            expandScalar(ff, f.asScalar());
+
+            forAll (ff, ffI)
+            {
+                localF[addr[ffI]] -= ff[ffI];
+            }
+        }
+    }
+    else if (this->activeType() == blockCoeffBase::SCALAR)
+    {
+        if (f.activeType() == blockCoeffBase::SQUARE)
+        {
+            squareTypeField& localF = this->asSquare();
+
+            const squareTypeField& ff = f.asSquare();
+
+            forAll (ff, ffI)
+            {
+                localF[addr[ffI]] -= ff[ffI];
+            }
+        }
+        if (f.activeType() == blockCoeffBase::LINEAR)
+        {
+            linearTypeField& localF = this->asLinear();
+
+            const linearTypeField& ff = f.asLinear();
+
+            forAll (ff, ffI)
+            {
+                localF[addr[ffI]] -= ff[ffI];
+            }
+        }
+        else if (f.activeType() == blockCoeffBase::SCALAR)
+        {
+            const scalarTypeField& ff = f.asScalar();
+
+            scalarTypeField& localF = this->asScalar();
+
+            forAll (ff, ffI)
+            {
+                localF[addr[ffI]] -= ff[ffI];
+            }
+        }
+    }
+    else
+    {
+        if (f.activeType() == blockCoeffBase::SQUARE)
+        {
+            squareTypeField& localF = this->toSquare();
+
+            const squareTypeField& ff = f.asSquare();
+
+            forAll (ff, ffI)
+            {
+                localF[addr[ffI]] -= ff[ffI];
+            }
+        }
+        if (f.activeType() == blockCoeffBase::LINEAR)
+        {
+            linearTypeField& localF = this->toLinear();
+
+            const linearTypeField& ff = f.asLinear();
+
+            forAll (ff, ffI)
+            {
+                localF[addr[ffI]] -= ff[ffI];
+            }
+        }
+        else if (f.activeType() == blockCoeffBase::SCALAR)
+        {
+            scalarTypeField& localF = this->toScalar();
+
+            const scalarTypeField& ff = f.asScalar();
+
+            forAll (ff, ffI)
+            {
+                localF[addr[ffI]] -= ff[ffI];
+            }
+        }
+//         FatalErrorIn
+//         (
+//             "template<class Type>\n"
+//             "void Foam::CoeffField<Type>::subtractSubset\n"
+//             "(\n"
+//             "    const CoeffField<Type>& f,\n"
+//             "    const labelList addr\n"
+//             ")"
+//         )   << "Incompatible combination of types"
+//             << abort(FatalError);
+    }
+}
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
@@ -1056,7 +1270,7 @@ void Foam::CoeffField<Type>::operator=(const CoeffField<Type>& f)
 
 
 template<class Type>
-void Foam::CoeffField<Type>::operator=(const tmp<CoeffField>& tf)
+void Foam::CoeffField<Type>::operator=(const tmp<CoeffField<Type> >& tf)
 {
     if (this == &(tf()))
     {
@@ -1070,7 +1284,16 @@ void Foam::CoeffField<Type>::operator=(const tmp<CoeffField>& tf)
 }
 
 
-#define COMPUTED_BASE_ASSIGNMENT(op)                                          \
+template<class Type>
+Foam::tmp<Foam::CoeffField<Type> > Foam::CoeffField<Type>::operator-()
+{
+    tmp<CoeffField<Type> > tf(new CoeffField<Type>(*this));
+    tf->negate();
+    return tf;
+}
+
+
+#define COMPUTED_BASE_ASSIGNMENT(op,defaultOp)                                \
                                                                               \
 template<class Type>                                                          \
 void Foam::CoeffField<Type>::operator op(const CoeffField<Type>& f)           \
@@ -1140,6 +1363,21 @@ void Foam::CoeffField<Type>::operator op(const CoeffField<Type>& f)           \
             this->asScalar() op f.asScalar();                                 \
         }                                                                     \
     }                                                                         \
+    else                                                                      \
+    {                                                                         \
+        if (f.activeType() == blockCoeffBase::SQUARE)                         \
+        {                                                                     \
+            (*this) defaultOp f.asSquare();                                   \
+        }                                                                     \
+        if (f.activeType() == blockCoeffBase::LINEAR)                         \
+        {                                                                     \
+            (*this) defaultOp f.asLinear();                                   \
+        }                                                                     \
+        else if (f.activeType() == blockCoeffBase::SCALAR)                    \
+        {                                                                     \
+            (*this) defaultOp f.asScalar();                                   \
+        }                                                                     \
+    }                                                                         \
 }                                                                             \
                                                                               \
 template<class Type>                                                          \
@@ -1173,9 +1411,6 @@ void Foam::CoeffField<Type>::operator op(const scalarTypeField& f)            \
         expandScalar(stf, f);                                                 \
         this->toSquare() op stf;                                              \
     }                                                                         \
-    else                                                                      \
-    {                                                                         \
-    }                                                                         \
 }                                                                             \
                                                                               \
 template<class Type>                                                          \
@@ -1208,9 +1443,6 @@ void Foam::CoeffField<Type>::operator op(const linearTypeField& f)            \
         expandLinear(stf, f);                                                 \
         this->toSquare() op stf;                                              \
     }                                                                         \
-    else                                                                      \
-    {                                                                         \
-    }                                                                         \
 }                                                                             \
                                                                               \
 template<class Type>                                                          \
@@ -1287,15 +1519,15 @@ void Foam::CoeffField<Type>::operator op(const tmp<Field<TYPE> >& tf)         \
 }
 
 
-#define COMPUTED_ASSIGNMENT(op)                                               \
+#define COMPUTED_ASSIGNMENT(op, defaultOp)                                    \
-COMPUTED_BASE_ASSIGNMENT(op)                                                  \
+COMPUTED_BASE_ASSIGNMENT(op, defaultOp)                                       \
 COMPUTED_ARG_ASSIGNMENT(op)
 
 // Remaining operator=
 COMPUTED_ARG_ASSIGNMENT(=)
 
-COMPUTED_ASSIGNMENT(+=)
+COMPUTED_ASSIGNMENT(+=, =)
-COMPUTED_ASSIGNMENT(-=)
+COMPUTED_ASSIGNMENT(-=, =-)
 
 COMPUTED_BASE_OPERATOR(scalar, *=)
 COMPUTED_BASE_OPERATOR(scalar, /=)

src/OpenFOAM/fields/CoeffField/CoeffField.H

@@ -252,7 +252,12 @@ public:
                 const labelList& addr
             );
 
-
+            //- Subtract subset with addressing to field
+            void subtractSubset
+            (
+                const CoeffField<Type>& f,
+                const labelList& addr
+            );
     // Member operators
 
         void operator=(const CoeffField<Type>&);
@@ -294,6 +299,7 @@ public:
         void operator/=(const tmp<Field<scalar> >&);
         void operator/=(const scalar&);
 
+        tmp<CoeffField<Type> > operator-();
 
     // IOstream operators
 

src/OpenFOAM/fields/CoeffField/CoeffFieldFunctions.C

@@ -327,9 +327,7 @@ Foam::tmp<Foam::Field<Type> > Foam::operator op                               \
     BINARY_OPERATOR_TRF(Type1, Type2, op, opFunc)                             \
     BINARY_OPERATOR_TRT(Type1, Type2, op, opFunc)
 
-// Operator multiply is not available for all types, as it expands rank
+BINARY_OPERATOR_R(Type, Type, *, multiply)
-// HJ, 17/Jun/2010
-// BINARY_OPERATOR_R(Type, Type, *, multiply)
 
 #undef BINARY_OPERATOR_R
 #undef BINARY_OPERATOR_FF

src/OpenFOAM/fields/CoeffField/DecoupledCoeffFieldFunctions.C

@@ -313,7 +313,7 @@ Foam::tmp<Foam::Field<Type> > Foam::operator op                               \
     BINARY_OPERATOR_TRF(Type1, Type2, op, opFunc)                             \
     BINARY_OPERATOR_TRT(Type1, Type2, op, opFunc)
 
-// BINARY_OPERATOR_R(Type, Type, *, multiply)
+BINARY_OPERATOR_R(Type, Type, *, multiply)
 
 #undef BINARY_OPERATOR_R
 #undef BINARY_OPERATOR_FF

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.C

@@ -397,6 +397,29 @@ interfaces() const
 }
 
 
+template<class Type, template<class> class PatchField, class GeoMesh>
+typename Foam::BlockLduInterfaceFieldPtrsList<Type>::Type
+Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricBoundaryField::
+blockInterfaces() const
+{
+    typename BlockLduInterfaceFieldPtrsList<Type>::Type interfaces(this->size());
+
+    forAll (interfaces, patchi)
+    {
+        if (isA<BlockLduInterfaceField<Type> >(this->operator[](patchi)))
+        {
+            interfaces.set
+            (
+                patchi,
+                &refCast<const BlockLduInterfaceField<Type> >(this->operator[](patchi))
+            );
+        }
+    }
+
+    return interfaces;
+}
+
+
 template<class Type, template<class> class PatchField, class GeoMesh>
 void Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricBoundaryField::
 writeEntry(const word& keyword, Ostream& os) const

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H

@@ -45,6 +45,7 @@ SourceFiles
 #include "DimensionedField.H"
 #include "FieldField.H"
 #include "lduInterfaceFieldPtrsList.H"
+#include "BlockLduInterfaceFieldPtrsList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -192,6 +193,10 @@ public:
             //  pointing to interfaces being set
             lduInterfaceFieldPtrsList interfaces() const;
 
+            //- Return a list of pointers for each patch field with only those
+            //  pointing to block-coupled interfaces being set
+            typename BlockLduInterfaceFieldPtrsList<Type>::Type blockInterfaces() const;
+            
             //- Write boundary field as dictionary entry
             void writeEntry(const word& keyword, Ostream& os) const;
  

src/OpenFOAM/matrices/blockLduMatrix/BlockLduMatrix/BlockConstraint/BlockConstraintTools.C

@@ -247,14 +247,14 @@ void BlockConstraint<Type>::setSourceDiag
     if (mag(fc) > SMALL)
     {
         b[rowID()] =
-            scale
+            cmptMultiply
             (
                 fc,
                 mult(matrix.diag().getCoeff(rowID()), value())
             );
 
         // set the solution to the right value as well
-        x[rowID()] = scale(fc, value());
+        x[rowID()] = cmptMultiply(fc, value());
     }
 }
 

src/OpenFOAM/matrices/blockLduMatrix/BlockLduMatrix/BlockLduInterface/BlockLduInterfaceField.H

@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+\    /   O peration     |
+ \  /    A nd           | Copyright held by original author
+  \/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    BlockLduInterfaceField
+
+Description
+    Abstract base class for interface fields with block coefficients
+
+Author
+    Ivor Clifford, ivor.clifford@gmail.com
+
+SourceFiles
+    BlockLduInterfaceField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef BlockLduInterfaceField_H
+#define BlockLduInterfaceField_H
+
+#include "lduInterfaceField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+class BlockLduMatrix;
+
+template<class Type>
+class CoeffField;
+
+/*---------------------------------------------------------------------------*\
+                           Class BlockLduInterfaceField Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class BlockLduInterfaceField
+:
+    public lduInterfaceField
+{
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        BlockLduInterfaceField(const BlockLduInterfaceField&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const BlockLduInterfaceField&);
+
+
+public:
+
+    //- Runtime Type information
+    TypeName("BlockLduInterfaceField");
+
+    // Constructors
+
+        //- Construct given coupled patch
+        BlockLduInterfaceField(const lduInterface& patch)
+        :
+            lduInterfaceField(patch)
+        {}
+
+
+    // Destructor
+
+        virtual ~BlockLduInterfaceField()
+        {}
+
+
+
+    // Member Functions
+
+        // Coupled interface matrix update
+
+            //- Initialise neighbour matrix update
+            virtual void initInterfaceMatrixUpdate
+            (
+                const Field<Type>&,
+                Field<Type>&,
+                const BlockLduMatrix<Type>&,
+                const CoeffField<Type>&,
+                const Pstream::commsTypes commsType
+            ) const
+            {}
+
+            //- Update result field based on interface functionality
+            virtual void updateInterfaceMatrix
+            (
+                const Field<Type>&,
+                Field<Type>&,
+                const BlockLduMatrix<Type>&,
+                const CoeffField<Type>&,
+                const Pstream::commsTypes commsType
+            ) const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/matrices/blockLduMatrix/BlockLduMatrix/BlockLduInterface/BlockLduInterfaceFieldPtrsList.H

@@ -0,0 +1,66 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    BlockLduInterfaceFieldPtrsList
+
+Description
+    List of coupled interface fields to be used in coupling.
+    Templated typedef workaround.
+
+Author
+    Ivor Clifford, ivor.clifford@gmail.com
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef BlockLduInterfaceFieldPtrsList_H
+#define BlockLduInterfaceFieldPtrsList_H
+
+#include "BlockLduInterfaceField.H"
+#include "UPtrList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<class T>
+struct BlockLduInterfaceFieldPtrsList
+{
+    typedef UPtrList<const BlockLduInterfaceField<T> > Type;
+};
+
+template<class T>
+struct BlockLduInterfaceFieldPtrsListList
+{
+    typedef List<typename BlockLduInterfaceFieldPtrsList<T>::Type> Type;
+};
+  
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/matrices/blockLduMatrix/BlockLduMatrix/BlockLduInterface/BlockLduInterfaceFields.C

@@ -0,0 +1,55 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    BlockLduInterfaceField
+
+Description
+    Abstract base class for interface fields with block coefficients
+
+Author
+    Ivor Clifford, ivor.clifford@gmail.com
+
+\*----------------------------------------------------------------------------*/
+
+#include "BlockLduInterfaceField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTemplateTypeNameAndDebug(BlockLduInterfaceField<scalar>, 0);
+defineTemplateTypeNameAndDebug(BlockLduInterfaceField<vector>, 0);
+defineTemplateTypeNameAndDebug(BlockLduInterfaceField<tensor>, 0);
+defineTemplateTypeNameAndDebug(BlockLduInterfaceField<sphericalTensor>, 0);
+defineTemplateTypeNameAndDebug(BlockLduInterfaceField<symmTensor>, 0);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/OpenFOAM/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrix.C

@@ -54,19 +54,8 @@ Foam::BlockLduMatrix<Type>::BlockLduMatrix(const lduMesh& ldu)
     diagPtr_(NULL),
     upperPtr_(NULL),
     lowerPtr_(NULL),
-    interfaces_(),
-    coupleUpper_(ldu.lduAddr().nPatches()),
-    coupleLower_(ldu.lduAddr().nPatches()),
     fixedEqns_(ldu.lduAddr().size()/fixFillIn)
-{
+{}
-    const lduAddressing& addr = ldu.lduAddr();
-
-    forAll (coupleUpper_, i)
-    {
-        coupleUpper_.set(i, new CoeffField<Type>(addr.patchAddr(i).size()));
-        coupleLower_.set(i, new CoeffField<Type>(addr.patchAddr(i).size()));
-    }
-}
 
 
 template<class Type>
@@ -77,9 +66,6 @@ Foam::BlockLduMatrix<Type>::BlockLduMatrix(const BlockLduMatrix<Type>& A)
     diagPtr_(NULL),
     upperPtr_(NULL),
     lowerPtr_(NULL),
-    interfaces_(),
-    coupleUpper_(),
-    coupleLower_(),
     fixedEqns_(A.fixedEqns_)
 {
     if (A.diagPtr_)
@@ -96,12 +82,56 @@ Foam::BlockLduMatrix<Type>::BlockLduMatrix(const BlockLduMatrix<Type>& A)
     {
         lowerPtr_ = new TypeCoeffField(*(A.lowerPtr_));
     }
-
-    // Interface data
-    coupleUpper_ = A.coupleUpper_;
-    coupleLower_ = A.coupleUpper_;
 }
 
+template<class Type>
+Foam::BlockLduMatrix<Type>::BlockLduMatrix(BlockLduMatrix<Type>& A, bool reUse)
+:
+    refCount(),
+    lduMesh_(A.lduMesh_),
+    diagPtr_(NULL),
+    upperPtr_(NULL),
+    lowerPtr_(NULL),
+    fixedEqns_(A.fixedEqns_)
+{
+    if (reUse)
+    {
+        if (A.lowerPtr_)
+        {
+            lowerPtr_ = A.lowerPtr_;
+            A.lowerPtr_ = NULL;
+        }
+
+        if (A.diagPtr_)
+        {
+            diagPtr_ = A.diagPtr_;
+            A.diagPtr_ = NULL;
+        }
+
+        if (A.upperPtr_)
+        {
+            upperPtr_ = A.upperPtr_;
+            A.upperPtr_ = NULL;
+        }
+    }
+    else
+    {
+        if (A.diagPtr_)
+        {
+            diagPtr_ = new TypeCoeffField(*(A.diagPtr_));
+        }
+
+        if (A.upperPtr_)
+        {
+            upperPtr_ = new TypeCoeffField(*(A.upperPtr_));
+        }
+
+        if (A.lowerPtr_)
+        {
+            lowerPtr_ = new TypeCoeffField(*(A.lowerPtr_));
+        }
+    }
+}
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
@@ -291,39 +321,37 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const BlockLduMatrix<Type>& ldum)
 {
     if (ldum.diagPtr_)
     {
-        os  << *ldum.diagPtr_ << nl;
+        os.writeKeyword("diagonal") << *ldum.diagPtr_ <<  token::END_STATEMENT << nl;
     }
     else
     {
         // Dummy write for consistency
-        os  << typename BlockLduMatrix<Type>::TypeCoeffField
+        os.writeKeyword("diagonal")  << typename BlockLduMatrix<Type>::TypeCoeffField
-            (ldum.lduAddr().size()) << nl;
+            (ldum.lduAddr().size()) << token::END_STATEMENT << nl;
     }
 
     if (ldum.upperPtr_)
     {
-        os  << *ldum.upperPtr_ << nl;
+        os.writeKeyword("upper")  << *ldum.upperPtr_ <<  token::END_STATEMENT << nl;
     }
     else
     {
         // Dummy write for consistency
-        os  << typename BlockLduMatrix<Type>::TypeCoeffField
+        os.writeKeyword("upper")  << typename BlockLduMatrix<Type>::TypeCoeffField
-            (ldum.lduAddr().lowerAddr().size()) << nl;
+            (ldum.lduAddr().lowerAddr().size()) <<  token::END_STATEMENT << nl;
     }
 
     if (ldum.lowerPtr_)
     {
-        os  << *ldum.lowerPtr_ << nl;
+        os.writeKeyword("lower")  << *ldum.lowerPtr_ <<  token::END_STATEMENT << nl;
     }
     else
     {
         // Dummy write for consistency
-        os  << typename BlockLduMatrix<Type>::TypeCoeffField
+        os.writeKeyword("lower")  << typename BlockLduMatrix<Type>::TypeCoeffField
-            (ldum.lduAddr().lowerAddr().size()) << nl;
+            (ldum.lduAddr().lowerAddr().size()) <<  token::END_STATEMENT << nl;
     }
 
-    os << endl;
-
     os.check("Ostream& operator<<(Ostream&, const BlockLduMatrix<Type>&");
 
     return os;

src/OpenFOAM/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrix.H

@@ -49,8 +49,9 @@ SourceFiles
 #define BlockLduMatrix_H
 
 #include "coeffFields.H"
+#include "FieldField.H"
 #include "lduMesh.H"
-#include "BlockLduInterface.H"
+#include "BlockLduInterfaceFieldPtrsList.H"
 #include "Map.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -87,7 +88,6 @@ public:
         typedef Field<Type> TypeField;
         typedef BlockConstraint<Type> ConstraintType;
 
-
 private:
 
     // Private data
@@ -96,7 +96,7 @@ private:
         const lduMesh& lduMesh_;
 
 
-        // Block matrix elements
+        // Block matrix elements (excluding interfaces)
 
             //- Diagonal coefficients
             CoeffField<Type>* diagPtr_;
@@ -108,18 +108,6 @@ private:
             CoeffField<Type> *lowerPtr_;
 
 
-        // Coupling
-
-            //- List of coupled interfaces
-            PtrList<BlockLduInterface<Type> > interfaces_;
-
-            //- Coupled interface coefficients, upper
-            FieldField<CoeffField, Type> coupleUpper_;
-
-            //- Coupled interface coefficients, lower
-            FieldField<CoeffField, Type> coupleLower_;
-
-
         // Constraints
 
             //- Equation triangle map
@@ -216,6 +204,8 @@ public:
         //- Construct as copy
         BlockLduMatrix(const BlockLduMatrix<Type>&);
 
+        //- Construct as copy or re-use as specified.
+        BlockLduMatrix(BlockLduMatrix<Type>&, bool reUse);
 
     // Destructor
 
@@ -267,30 +257,6 @@ public:
             //- Return lower coefficients
             const TypeCoeffField& lower() const;
 
-            //- Return access to coupled interface coefficients, upper
-            FieldField<CoeffField, Type>& coupleUpper()
-            {
-                return coupleUpper_;
-            }
-
-            //- Return coupled interface coefficients, upper
-            const FieldField<CoeffField, Type>& coupleUpper() const
-            {
-                return coupleUpper_;
-            }
-
-            //- Return access to coupled interface coefficients, lower
-            FieldField<CoeffField, Type>& coupleLower()
-            {
-                return coupleLower_;
-            }
-
-            //- Return coupled interface coefficients, lower
-            const FieldField<CoeffField, Type>& coupleLower() const
-            {
-                return coupleLower_;
-            }
-
 
         // Matrix structure
 
@@ -351,7 +317,9 @@ public:
             void Amul
             (
                 TypeField& Ax,
-                const TypeField& x
+                const TypeField& x,
+				const FieldField<CoeffField, Type>& boundaryCoeffs,
+                const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces
             ) const;
 
             //- Matrix multiplication without coupled interface update
@@ -365,11 +333,12 @@ public:
             void Tmul
             (
                 TypeField& Ax,
-                const TypeField& x
+                const TypeField& x,
+				const FieldField<CoeffField, Type>& boundaryCoeffs,
+                const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces
             ) const;
 
-            //- Matrix transpose multiplication without
+            //- Matrix transpose multiplication without coupled interface update
-            //  coupled interface update
             void TmulCore
             (
                 TypeField& Ax,
@@ -391,16 +360,19 @@ public:
             //  for Amul operations
             void initInterfaces
             (
-                TypeField& Ax,
+                const FieldField<CoeffField, Type>& interfaceCoeffs,
-                const TypeField& x
+                const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
+                TypeField& result,
+                const TypeField& psi
             ) const;
 
             //- Update coupled interfaces
             void updateInterfaces
             (
-                const FieldField<CoeffField, Type>& coeffs,
+                const FieldField<CoeffField, Type>& interfaceCoeffs,
-                TypeField& Ax,
+                const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
-                const TypeField& x
+                TypeField& result,
+                const TypeField& psi
             ) const;
 
 

src/OpenFOAM/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixATmul.C

@@ -35,16 +35,20 @@ template<class Type>
 void Foam::BlockLduMatrix<Type>::Amul
 (
     TypeField& Ax,
-    const TypeField& x
+    const TypeField& x,
+    const FieldField<CoeffField, Type>& boundaryCoeffs,
+    const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces
 ) const
 {
+    Ax = pTraits<Type>::zero;
+
     // Initialise the update of coupled interfaces
-    initInterfaces(Ax, x);
+    initInterfaces(boundaryCoeffs, interfaces, Ax, x);
 
     AmulCore(Ax, x);
 
     // Update coupled interfaces
-    updateInterfaces(coupleUpper_, Ax, x);
+    updateInterfaces(boundaryCoeffs, interfaces, Ax, x);
 }
 
 
@@ -175,16 +179,20 @@ template<class Type>
 void Foam::BlockLduMatrix<Type>::Tmul
 (
     TypeField& Ax,
-    const TypeField& x
+    const TypeField& x,
+    const FieldField<CoeffField, Type>& boundaryCoeffs,
+    const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces
 ) const
 {
+    Ax = pTraits<Type>::zero;
+
     // Initialise the update of coupled interfaces
-    initInterfaces(Ax, x);
+    initInterfaces(boundaryCoeffs, interfaces, Ax, x);
 
     TmulCore(Ax, x);
 
     // Update coupled interfaces
-    updateInterfaces(coupleLower_, Ax, x);
+    updateInterfaces(boundaryCoeffs, interfaces, Ax, x);
 }
 
 

src/OpenFOAM/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixHOps.C

@@ -0,0 +1,273 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-6 H. Jasak All rights reserved
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+Description
+    Vector-matrix multiplication operations for a block matrix
+
+\*---------------------------------------------------------------------------*/
+
+#include "BlockLduMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::BlockLduMatrix<Type>::H(const Field<Type>& x) const
+{
+    typedef typename TypeCoeffField::scalarTypeField scalarTypeField;
+    typedef typename TypeCoeffField::linearTypeField linearTypeField;
+    typedef typename TypeCoeffField::squareTypeField squareTypeField;
+
+    // Create result
+    tmp<Field<Type> > tresult
+    (
+        new Field<Type>(lduAddr().size(), pTraits<Type>::zero)
+    );
+    Field<Type>& result = tresult();
+
+    if (this->thereIsUpper())
+    {
+        const unallocLabelList& l = lduAddr().lowerAddr();
+        const unallocLabelList& u = lduAddr().upperAddr();
+        const TypeCoeffField& Upper = this->upper();
+
+        // Create multiplication function object
+        typename BlockCoeff<Type>::multiply mult;
+
+        // Lower multiplication
+
+        if (symmetric())
+        {
+            if (Upper.activeType() == blockCoeffBase::SCALAR)
+            {
+                const scalarTypeField& activeUpper = Upper.asScalar();
+
+                for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+                {
+                    result[u[coeffI]] -= mult(activeUpper[coeffI], x[l[coeffI]]);
+                }
+            }
+            else if (Upper.activeType() == blockCoeffBase::LINEAR)
+            {
+                const linearTypeField& activeUpper = Upper.asLinear();
+
+                for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+                {
+                    result[u[coeffI]] -= mult(activeUpper[coeffI], x[l[coeffI]]);
+                }
+            }
+            else if (Upper.activeType() == blockCoeffBase::SQUARE)
+            {
+                const squareTypeField& activeUpper = Upper.asSquare();
+
+                for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+                {
+                    // Use transpose upper coefficient
+                    result[u[coeffI]] -=
+                        mult(activeUpper[coeffI].T(), x[l[coeffI]]);
+                }
+            }
+        }
+        else // Asymmetric matrix
+        {
+            const TypeCoeffField& Lower = this->lower();
+
+            if (Lower.activeType() == blockCoeffBase::SCALAR)
+            {
+                const scalarTypeField& activeLower = Lower.asScalar();
+
+                for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+                {
+                    result[u[coeffI]] -= mult(activeLower[coeffI], x[l[coeffI]]);
+                }
+            }
+            else if (Lower.activeType() == blockCoeffBase::LINEAR)
+            {
+                const linearTypeField& activeLower = Lower.asLinear();
+
+                for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+                {
+                    result[u[coeffI]] -= mult(activeLower[coeffI], x[l[coeffI]]);
+                }
+            }
+            else if (Lower.activeType() == blockCoeffBase::SQUARE)
+            {
+                const squareTypeField& activeLower = Lower.asSquare();
+
+                for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+                {
+                    result[u[coeffI]] -= mult(activeLower[coeffI], x[l[coeffI]]);
+                }
+            }
+        }
+
+
+        // Upper multiplication
+
+        if (Upper.activeType() == blockCoeffBase::SCALAR)
+        {
+            const scalarTypeField& activeUpper = Upper.asScalar();
+
+            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            {
+                result[l[coeffI]] -= mult(activeUpper[coeffI], x[u[coeffI]]);
+            }
+        }
+        else if (Upper.activeType() == blockCoeffBase::LINEAR)
+        {
+            const linearTypeField& activeUpper = Upper.asLinear();
+
+            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            {
+                result[l[coeffI]] -= mult(activeUpper[coeffI], x[u[coeffI]]);
+            }
+        }
+        else if (Upper.activeType() == blockCoeffBase::SQUARE)
+        {
+            const squareTypeField& activeUpper = Upper.asSquare();
+
+            for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+            {
+                result[l[coeffI]] -= mult(activeUpper[coeffI], x[u[coeffI]]);
+            }
+        }
+    }
+
+    return tresult;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::BlockLduMatrix<Type>::faceH(const Field<Type>& x) const
+{
+    typedef typename TypeCoeffField::scalarTypeField scalarTypeField;
+    typedef typename TypeCoeffField::linearTypeField linearTypeField;
+    typedef typename TypeCoeffField::squareTypeField squareTypeField;
+
+    const unallocLabelList& u = lduAddr().upperAddr();
+    const unallocLabelList& l = lduAddr().lowerAddr();
+
+    // Create result
+    tmp<Field<Type> > tresult(new Field<Type>(u.size(), pTraits<Type>::zero));
+    Field<Type>& result = tresult();
+
+    if (this->thereIsUpper())
+    {
+        const TypeCoeffField& Upper = this->upper();
+
+        // Create multiplication function object
+        typename BlockCoeff<Type>::multiply mult;
+
+        // Lower multiplication
+
+        if (symmetric())
+        {
+            if (Upper.activeType() == blockCoeffBase::SCALAR)
+            {
+                const scalarTypeField& activeUpper = Upper.asScalar();
+
+                for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+                {
+                    // This can be optimised with a subtraction.
+                    // Currently not done for clarity.  HJ, 31/Oct/2007
+                    result[coeffI] =
+                        mult(activeUpper[coeffI], x[u[coeffI]])
+                    - mult(activeUpper[coeffI], x[l[coeffI]]);
+                }
+            }
+            else if (Upper.activeType() == blockCoeffBase::LINEAR)
+            {
+                const linearTypeField& activeUpper = Upper.asLinear();
+
+                for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+                {
+                    // This can be optimised with a subtraction.
+                    // Currently not done for clarity.  HJ, 31/Oct/2007
+                    result[coeffI] =
+                        mult(activeUpper[coeffI], x[u[coeffI]])
+                    - mult(activeUpper[coeffI], x[l[coeffI]]);
+                }
+            }
+            else if (Upper.activeType() == blockCoeffBase::SQUARE)
+            {
+                const squareTypeField& activeUpper = Upper.asSquare();
+
+                for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+                {
+                    // Use transpose upper coefficient
+                    result[coeffI] =
+                        mult(activeUpper[coeffI], x[u[coeffI]])
+                    - mult(activeUpper[coeffI].T(), x[l[coeffI]]);
+                }
+            }
+        }
+        else // Asymmetric matrix
+        {
+            const TypeCoeffField& Lower = this->lower();
+
+            if (Lower.activeType() == blockCoeffBase::SCALAR)
+            {
+                const scalarTypeField& activeUpper = Upper.asScalar();
+                const scalarTypeField& activeLower = Lower.asScalar();
+
+                for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+                {
+                    result[coeffI] =
+                        mult(activeUpper[coeffI], x[u[coeffI]])
+                    - mult(activeLower[coeffI], x[l[coeffI]]);
+                }
+            }
+            else if (Lower.activeType() == blockCoeffBase::LINEAR)
+            {
+                const linearTypeField& activeUpper = Upper.asLinear();
+                const linearTypeField& activeLower = Lower.asLinear();
+
+                for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+                {
+                    result[coeffI] =
+                        mult(activeUpper[coeffI], x[u[coeffI]])
+                    - mult(activeLower[coeffI], x[l[coeffI]]);
+                }
+            }
+            else if (Lower.activeType() == blockCoeffBase::SQUARE)
+            {
+                const squareTypeField& activeUpper = Upper.asSquare();
+                const squareTypeField& activeLower = Lower.asSquare();
+
+                for (register label coeffI = 0; coeffI < u.size(); coeffI++)
+                {
+                    result[coeffI] =
+                        mult(activeUpper[coeffI], x[u[coeffI]])
+                    - mult(activeLower[coeffI], x[l[coeffI]]);
+                }
+            }
+        }
+    }
+
+    return tresult;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixOperations.C

@@ -805,10 +805,6 @@ void Foam::BlockLduMatrix<Type>::negate()
     {
         diagPtr_->negate();
     }
-
-    // Do coupling coefficients
-    coupleUpper_.negate();
-    coupleLower_.negate();
 }
 
 
@@ -853,11 +849,6 @@ void Foam::BlockLduMatrix<Type>::operator=(const BlockLduMatrix<Type>& A)
         diag() = A.diag();
     }
 
-    // Copy interface data
-    interfaces_ = A.interfaces_;
-    coupleUpper_ = A.coupleUpper_;
-    coupleLower_ = A.coupleLower_;
-
     // Copy constraints
     fixedEqns_ = A.fixedEqns_;
 }
@@ -887,10 +878,6 @@ void Foam::BlockLduMatrix<Type>::operator+=(const BlockLduMatrix<Type>& A)
         upper() += A.upper();
         lower() += A.lower();
     }
-
-    // Interface data
-    coupleUpper_ += A.coupleUpper_;
-    coupleLower_ += A.coupleLower_;
 }
 
 
@@ -918,10 +905,6 @@ void Foam::BlockLduMatrix<Type>::operator-=(const BlockLduMatrix<Type>& A)
         upper() -= A.upper();
         lower() -= A.lower();
     }
-
-    // Interface data
-    coupleUpper_ -= A.coupleUpper_;
-    coupleLower_ -= A.coupleLower_;
 }
 
 
@@ -934,6 +917,9 @@ void Foam::BlockLduMatrix<Type>::operator*=(const scalarField& sf)
 
     //HJ Missing code: add coupling coefficients op*=
     //   HJ, 21/Feb/2008
+    // IC - Complicated because we have to scale across the interfaces
+    // We may need to include this functionality in lduInterfaceField
+    // as additional initInterfaceScale and scaleInterface functions
 
     if (diagPtr_)
     {
@@ -1029,10 +1015,6 @@ void Foam::BlockLduMatrix<Type>::operator*=(const scalar s)
     {
         *lowerPtr_ *= s;
     }
-
-    // Interface data
-    coupleUpper_ *= s;
-    coupleLower_ *= s;
 }
 
 

src/OpenFOAM/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixUpdateInterfaces.C

@@ -0,0 +1,204 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-6 H. Jasak All rights reserved
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+Description
+    Update of block interfaces
+
+\*---------------------------------------------------------------------------*/
+
+#include "BlockLduMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::BlockLduMatrix<Type>::initInterfaces
+(
+    const FieldField<CoeffField, Type>& interfaceCoeffs,
+    const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
+    TypeField& result,
+    const TypeField& psi
+) const
+{
+    if
+    (
+        Pstream::defaultCommsType == Pstream::blocking
+     || Pstream::defaultCommsType == Pstream::nonBlocking
+    )
+    {
+        forAll (interfaces, interfaceI)
+        {
+            if (interfaces.set(interfaceI))
+            {
+                interfaces[interfaceI].initInterfaceMatrixUpdate
+                (
+                    psi,
+                    result,
+                    *this,
+                    interfaceCoeffs[interfaceI],
+                    Pstream::defaultCommsType
+                );
+            }
+        }
+    }
+    else if (Pstream::defaultCommsType == Pstream::scheduled)
+    {
+        const lduSchedule& patchSchedule = this->patchSchedule();
+
+        // Loop over the "global" patches are on the list of interfaces but
+        // beyond the end of the schedule which only handles "normal" patches
+        for
+        (
+            label interfaceI=patchSchedule.size()/2;
+            interfaceI<interfaces.size();
+            interfaceI++
+        )
+        {
+            if (interfaces.set(interfaceI))
+            {
+                interfaces[interfaceI].initInterfaceMatrixUpdate
+                (
+                    psi,
+                    result,
+                    *this,
+                    interfaceCoeffs[interfaceI],
+                    Pstream::blocking
+                );
+            }
+        }
+    }    
+    else
+    {
+        FatalErrorIn("BlockLduMatrix<Type>::initMatrixInterfaces")
+            << "Unsuported communications type "
+            << Pstream::commsTypeNames[Pstream::defaultCommsType]
+            << exit(FatalError);
+    }
+}
+
+
+template<class Type>
+void Foam::BlockLduMatrix<Type>::updateInterfaces
+(
+    const FieldField<CoeffField, Type>& interfaceCoeffs,
+    const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
+    TypeField& result,
+    const TypeField& psi
+) const
+{
+    if
+    (
+        Pstream::defaultCommsType == Pstream::blocking
+     || Pstream::defaultCommsType == Pstream::nonBlocking
+    )
+    {
+        // Block until all sends/receives have been finished
+        if (Pstream::defaultCommsType == Pstream::nonBlocking)
+        {
+            IPstream::waitRequests();
+            OPstream::waitRequests();
+        }
+
+        forAll (interfaces, interfaceI)
+        {
+            if (interfaces.set(interfaceI))
+            {
+                interfaces[interfaceI].updateInterfaceMatrix
+                (
+                    psi,
+                    result,
+                    *this,
+                    interfaceCoeffs[interfaceI],
+                    Pstream::defaultCommsType
+                );
+            }
+        }
+    }
+    else if (Pstream::defaultCommsType == Pstream::scheduled)
+    {
+        const lduSchedule& patchSchedule = this->patchSchedule();
+
+        // Loop over all the "normal" interfaces relating to standard patches
+        forAll (patchSchedule, i)
+        {
+            label interfaceI = patchSchedule[i].patch;
+
+            if (interfaces.set(interfaceI))
+            {
+                if (patchSchedule[i].init)
+                {
+                    interfaces[interfaceI].initInterfaceMatrixUpdate
+                    (
+                        psi,
+                        result,
+                        *this,
+                        interfaceCoeffs[interfaceI],
+                        Pstream::scheduled
+                    );
+                }
+                else
+                {
+                    interfaces[interfaceI].updateInterfaceMatrix
+                    (
+                        psi,
+                        result,
+                        *this,
+                        interfaceCoeffs[interfaceI],
+                        Pstream::scheduled
+                    );
+                }
+            }
+        }
+
+        // Loop over the "global" patches are on the list of interfaces but
+        // beyond the end of the schedule which only handles "normal" patches
+        for
+        (
+            label interfaceI=patchSchedule.size()/2;
+            interfaceI<interfaces.size();
+            interfaceI++
+        )
+        {
+            if (interfaces.set(interfaceI))
+            {
+                interfaces[interfaceI].updateInterfaceMatrix
+                (
+                    psi,
+                    result,
+                    *this,
+                    interfaceCoeffs[interfaceI],
+                    Pstream::blocking
+                );
+            }
+        }
+    }
+    else
+    {
+        FatalErrorIn("BlockLduMatrix<Type>::updateMatrixInterfaces")
+            << "Unsuported communications type "
+            << Pstream::commsTypeNames[Pstream::defaultCommsType]
+            << exit(FatalError);
+    }
+}
+
+// ************************************************************************* //

src/OpenFOAM/matrices/blockLduMatrix/BlockLduMatrix/scalarBlockLduMatrix.C

@@ -24,9 +24,6 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef scalarBlockLduMatrix_H
-#define scalarBlockLduMatrix_H
-
 #include "BlockLduMatrix.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -423,8 +420,4 @@ tmp<scalarField> BlockLduMatrix<scalar>::faceH(const scalarField& x) const
 
 } // End namespace Foam
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
 // ************************************************************************* //

src/OpenFOAM/matrices/blockLduMatrix/BlockLduMatrix/scalarBlockLduMatrixUpdateInterfaces.C

@@ -0,0 +1,206 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-6 H. Jasak All rights reserved
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+Description
+    Update of block interfaces
+
+\*---------------------------------------------------------------------------*/
+
+#include "BlockLduMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// TODO - This code is currently not called so we have specialized
+// initInterfaceMatrixUpdate in processorFvPatchScalarfield. This needs to be fixed
+
+template<>
+void Foam::BlockLduMatrix<scalar>::initInterfaces
+(
+    const FieldField<CoeffField, scalar>& coupleCoeffs,
+    scalarField& result,
+    const scalarField& psi
+) const
+{
+    if
+    (
+        Pstream::defaultCommsType == Pstream::blocking
+     || Pstream::defaultCommsType == Pstream::nonBlocking
+    )
+    {
+        forAll (interfaces, interfaceI)
+        {
+            if (interfaces.set(interfaceI))
+            {
+                interfaces[interfaceI].initInterfaceMatrixUpdate
+                (
+                    psi,
+                    result,
+                    *this,
+                    coupleCoeffs[interfaceI].asScalar(),
+                    Pstream::defaultCommsType
+                );
+            }
+        }
+    }
+    else if (Pstream::defaultCommsType == Pstream::scheduled)
+    {
+        const lduSchedule& patchSchedule = this->patchSchedule();
+
+        // Loop over the "global" patches are on the list of interfaces but
+        // beyond the end of the schedule which only handles "normal" patches
+        for
+        (
+            label interfaceI=patchSchedule.size()/2;
+            interfaceI<interfaces.size();
+            interfaceI++
+        )
+        {
+            if (interfaces.set(interfaceI))
+            {
+                interfaces[interfaceI].initInterfaceMatrixUpdate
+                (
+                    psi,
+                    result,
+                    *this,
+                    coupleCoeffs[interfaceI].asScalar(),
+                    Pstream::blocking
+                );
+            }
+        }
+    }    
+    else
+    {
+        FatalErrorIn("BlockLduMatrix<scalar>::initMatrixInterfaces")
+            << "Unsuported communications type "
+            << Pstream::commsTypeNames[Pstream::defaultCommsType]
+            << exit(FatalError);
+    }
+}
+
+
+template<>
+void Foam::BlockLduMatrix<scalar>::updateInterfaces
+(
+    const FieldField<CoeffField, scalar>& coupleCoeffs,
+    scalarField& result,
+    const scalarField& psi
+) const
+{
+    if
+    (
+        Pstream::defaultCommsType == Pstream::blocking
+     || Pstream::defaultCommsType == Pstream::nonBlocking
+    )
+    {
+        // Block until all sends/receives have been finished
+        if (Pstream::defaultCommsType == Pstream::nonBlocking)
+        {
+            IPstream::waitRequests();
+            OPstream::waitRequests();
+        }
+
+        forAll (interfaces, interfaceI)
+        {
+            if (interfaces.set(interfaceI))
+            {
+                interfaces[interfaceI].updateInterfaceMatrix
+                (
+                    psi,
+                    result,
+                    *this,
+                    coupleCoeffs[interfaceI].asScalar(),
+                    Pstream::defaultCommsType
+                );
+            }
+        }
+    }
+    else if (Pstream::defaultCommsType == Pstream::scheduled)
+    {
+        const lduSchedule& patchSchedule = this->patchSchedule();
+
+        // Loop over all the "normal" interfaces relating to standard patches
+        forAll (patchSchedule, i)
+        {
+            label interfaceI = patchSchedule[i].patch;
+
+            if (interfaces.set(interfaceI))
+            {
+                if (patchSchedule[i].init)
+                {
+                    interfaces[interfaceI].initInterfaceMatrixUpdate
+                    (
+                        psi,
+                        result,
+                        *this,
+                        coupleCoeffs[interfaceI].asScalar(),
+                        Pstream::scheduled
+                    );
+                }
+                else
+                {
+                    interfaces[interfaceI].updateInterfaceMatrix
+                    (
+                        psi,
+                        result,
+                        *this,
+                        coupleCoeffs[interfaceI].asScalar(),
+                        Pstream::scheduled
+                    );
+                }
+            }
+        }
+
+        // Loop over the "global" patches are on the list of interfaces but
+        // beyond the end of the schedule which only handles "normal" patches
+        for
+        (
+            label interfaceI=patchSchedule.size()/2;
+            interfaceI<interfaces.size();
+            interfaceI++
+        )
+        {
+            if (interfaces.set(interfaceI))
+            {
+                interfaces[interfaceI].updateInterfaceMatrix
+                (
+                    psi,
+                    result,
+                    *this,
+                    coupleCoeffs[interfaceI].asScalar(),
+                    Pstream::blocking
+                );
+            }
+        }
+    }
+    else
+    {
+        FatalErrorIn("BlockLduMatrix<scalar>::updateMatrixInterfaces")
+            << "Unsuported communications type "
+            << Pstream::commsTypeNames[Pstream::defaultCommsType]
+            << exit(FatalError);
+    }
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/blockLduMatrix/BlockLduPrecons/BlockLduPrecon/BlockLduPrecon.H

@@ -70,12 +70,6 @@ protected:
         const BlockLduMatrix<Type>& matrix_;
 
 
-    // Protected member functions
-
-        //- Find the smoother name (directly or from a sub-dictionary)
-        static word getName(const dictionary&);
-
-
 public:
 
     //- Runtime type information
@@ -157,7 +151,7 @@ public:
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-#   include "BlockLduPrecon.C"
+#   include "newBlockLduPrecon.C"
 #endif
 
 #define makeBlockPrecon(PreconType, typePreconType)                          \

src/OpenFOAM/matrices/blockLduMatrix/BlockLduPrecons/BlockLduPrecon/newBlockLduPrecon.C

@@ -22,7 +22,7 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-\*---------------------------------------------------------------------------*/
+\*----------------------------------------------------------------------------*/
 
 #include "BlockLduPrecon.H"
 #include "blockNoPrecons.H"
@@ -39,20 +39,7 @@ Foam::autoPtr<Foam::BlockLduPrecon<Type> > Foam::BlockLduPrecon<Type>::New
     const dictionary& dict
 )
 {
-    word preconName;
+    word preconName(dict.lookup("preconditioner"));
-
-    // handle primitive or dictionary entry
-    const entry& e = dict.lookupEntry("preconditioner", false, false);
-    if (e.isDict())
-    {
-        e.dict().lookup("preconditioner") >> preconName;
-    }
-    else
-    {
-        e.stream() >> preconName;
-    }
-
-    const dictionary& controls = e.isDict() ? e.dict() : dictionary::null;
 
     if (matrix.diagonal())
     {
@@ -100,26 +87,4 @@ Foam::autoPtr<Foam::BlockLduPrecon<Type> > Foam::BlockLduPrecon<Type>::New
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class Type>
-Foam::word Foam::BlockLduPrecon<Type>::getName(const dictionary& dict)
-{
-    word name;
-
-    // handle primitive or dictionary entry
-    const entry& e = dict.lookupEntry("preconditioner", false, false);
-    if (e.isDict())
-    {
-        e.dict().lookup("preconditioner") >> name;
-    }
-    else
-    {
-        e.stream() >> name;
-    }
-
-    return name;
-}
-
-
 // ************************************************************************* //

src/OpenFOAM/matrices/blockLduMatrix/BlockLduSmoothers/BlockGaussSeidelSmoother/BlockGaussSeidelSmoother.H

@@ -83,10 +83,12 @@ public:
         BlockGaussSeidelSmoother
         (
             const BlockLduMatrix<Type>& matrix,
+            const FieldField<CoeffField, Type>& boundaryCoeffs,
+            const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
             const dictionary& dict
         )
         :
-            BlockLduSmoother<Type>(matrix),
+            BlockLduSmoother<Type>(matrix, boundaryCoeffs, interfaces),
             gs_(matrix)
         {}
 

src/OpenFOAM/matrices/blockLduMatrix/BlockLduSmoothers/BlockILUSmoother/BlockILUSmoother.H

@@ -86,10 +86,12 @@ public:
         BlockILUSmoother
         (
             const BlockLduMatrix<Type>& matrix,
+            const FieldField<CoeffField, Type>& boundaryCoeffs,
+            const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
             const dictionary& dict
         )
         :
-            BlockLduSmoother<Type>(matrix),
+            BlockLduSmoother<Type>(matrix, boundaryCoeffs, interfaces),
             precon_(matrix),
             xCorr_(matrix.lduAddr().size()),
             residual_(matrix.lduAddr().size())
@@ -115,7 +117,13 @@ public:
             for (label sweep = 0; sweep < nSweeps; sweep++)
             {
                 // Calculate residual
-                this-> matrix_.Amul(residual_, x);
+                this-> matrix_.Amul
+                (
+                    residual_,
+                    x,
+                    BlockLduSmoother<Type>::boundaryCoeffs_,
+                    BlockLduSmoother<Type>::interfaces_
+                );
 
                 // residual = b - Ax
                 forAll (b, i)

src/OpenFOAM/matrices/blockLduMatrix/BlockLduSmoothers/BlockLduSmoother/BlockLduSmoother.H

@@ -66,12 +66,11 @@ protected:
         //- Matrix reference
         const BlockLduMatrix<Type>& matrix_;
 
+        //- Reference to boundary coefficients
+        const FieldField<CoeffField, Type>& boundaryCoeffs_;
 
-    // Protected member functions
+        //- Reference to interfaces
-
+        const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces_;
-        //- Find the smoother name (directly or from a sub-dictionary)
-        static word getName(const dictionary&);
-
 
 public:
 
@@ -88,10 +87,14 @@ public:
             dictionary,
             (
                 const BlockLduMatrix<Type>& matrix,
+                const FieldField<CoeffField, Type>& boundaryCoeffs,
+                const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
                 const dictionary& dict
             ),
             (
                 matrix,
+                boundaryCoeffs,
+                interfaces,
                 dict
             )
         );
@@ -100,9 +103,16 @@ public:
     // Constructors
 
         //- Construct from matrix
-        BlockLduSmoother(const BlockLduMatrix<Type>& matrix)
+        BlockLduSmoother
+        (
+            const BlockLduMatrix<Type>& matrix,
+            const FieldField<CoeffField, Type>& boundaryCoeffs,
+            const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces
+        )
         :
-            matrix_(matrix)
+            matrix_(matrix),
+            boundaryCoeffs_(boundaryCoeffs),
+            interfaces_(interfaces)
         {}
 
 
@@ -112,6 +122,8 @@ public:
         static autoPtr<BlockLduSmoother<Type> > New
         (
             const BlockLduMatrix<Type>& matrix,
+            const FieldField<CoeffField, Type>& boundaryCoeffs,
+            const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
             const dictionary& dict
         );
 
@@ -141,7 +153,7 @@ public:
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-#   include "BlockLduSmoother.C"
+#   include "newBlockLduSmoother.C"
 #endif
 
 #define makeBlockSmoother(SmootherType, typeSmootherType)                    \
@@ -153,8 +165,8 @@ addToRunTimeSelectionTable(SmootherType, typeSmootherType, dictionary);
 #define makeBlockSmoothers(smootherType)                                       \
                                                                              \
 makeBlockSmoother(blockScalarSmoother, smootherType##Scalar);                  \
-makeBlockSmoother(blockVectorSmoother, smootherType##Vector);                  \
+makeBlockSmoother(blockVectorSmoother, smootherType##Vector);
-makeBlockSmoother(blockTensorSmoother, smootherType##Tensor);
+// makeBlockSmoother(blockTensorSmoother, smootherType##Tensor);
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/OpenFOAM/matrices/blockLduMatrix/BlockLduSmoothers/BlockLduSmoother/blockLduSmoothers.C

@@ -38,11 +38,11 @@ namespace Foam
 
 defineNamedTemplateTypeNameAndDebug(blockScalarSmoother, 0);
 defineNamedTemplateTypeNameAndDebug(blockVectorSmoother, 0);
-defineNamedTemplateTypeNameAndDebug(blockTensorSmoother, 0);
+// defineNamedTemplateTypeNameAndDebug(blockTensorSmoother, 0);
 
 defineTemplateRunTimeSelectionTable(blockScalarSmoother, dictionary);
 defineTemplateRunTimeSelectionTable(blockVectorSmoother, dictionary);
-defineTemplateRunTimeSelectionTable(blockTensorSmoother, dictionary);
+// defineTemplateRunTimeSelectionTable(blockTensorSmoother, dictionary);
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/OpenFOAM/matrices/blockLduMatrix/BlockLduSmoothers/BlockLduSmoother/blockLduSmoothers.H

@@ -50,7 +50,7 @@ namespace Foam
 
 typedef BlockLduSmoother<scalar> blockScalarSmoother;
 typedef BlockLduSmoother<vector> blockVectorSmoother;
-typedef BlockLduSmoother<tensor> blockTensorSmoother;
+// typedef BlockLduSmoother<tensor> blockTensorSmoother;
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/OpenFOAM/matrices/blockLduMatrix/BlockLduSmoothers/BlockLduSmoother/newBlockLduSmoother.C

@@ -28,7 +28,7 @@ Class
 Description
     Block LDU matrix smoother virtual base class
 
-\*---------------------------------------------------------------------------*/
+\*----------------------------------------------------------------------------*/
 
 #include "BlockLduSmoother.H"
 
@@ -41,24 +41,12 @@ template<class Type>
 Foam::autoPtr<Foam::BlockLduSmoother<Type> > Foam::BlockLduSmoother<Type>::New
 (
     const BlockLduMatrix<Type>& matrix,
+    const FieldField<CoeffField, Type>& boundaryCoeffs,
+    const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
     const dictionary& dict
 )
 {
-    word smootherName;
+    word smootherName(dict.lookup("smoother"));
-
-    // Handle primitive or dictionary entry
-    const entry& e = dict.lookupEntry("smoother", false, false);
-    if (e.isDict())
-    {
-        e.dict().lookup("smoother") >> smootherName;
-    }
-    else
-    {
-        e.stream() >> smootherName;
-    }
-
-    // Not (yet?) needed:
-    // const dictionary& controls = e.isDict() ? e.dict() : dictionary::null;
 
     typename dictionaryConstructorTable::iterator constructorIter =
         dictionaryConstructorTablePtr_->find(smootherName);
@@ -70,6 +58,8 @@ Foam::autoPtr<Foam::BlockLduSmoother<Type> > Foam::BlockLduSmoother<Type>::New
             "autoPtr<BlockLduSmoother> BlockLduSmoother::New\n"
             "(\n"
             "    const BlockLduMatrix<Type>& matrix,\n"
+            "    const FieldField<CoeffField, Type>& boundaryCoeffs,\n"
+            "    const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,\n"
             "    const dictionary& dict\n"
             ")",
             dict
@@ -85,32 +75,12 @@ Foam::autoPtr<Foam::BlockLduSmoother<Type> > Foam::BlockLduSmoother<Type>::New
         constructorIter()
         (
             matrix,
+            boundaryCoeffs,
+            interfaces,
             dict
         )
     );
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class Type>
-Foam::word Foam::BlockLduSmoother<Type>::getName(const dictionary& dict)
-{
-    word name;
-
-    // handle primitive or dictionary entry
-    const entry& e = dict.lookupEntry("preconditioner", false, false);
-    if (e.isDict())
-    {
-        e.dict().lookup("preconditioner") >> name;
-    }
-    else
-    {
-        e.stream() >> name;
-    }
-
-    return name;
-}
-
-
 // ************************************************************************* //

src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockBiCGStab/BlockBiCGStabSolver.C

@@ -37,14 +37,18 @@ Foam::BlockBiCGStabSolver<Type>::BlockBiCGStabSolver
 (
     const word& fieldName,
     const BlockLduMatrix<Type>& matrix,
-    const dictionary& dict
+    const FieldField<CoeffField, Type>& boundaryCoeffs,
+    const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
+    const dictionary& solverDict
 )
 :
     BlockIterativeSolver<Type>
     (
         fieldName,
         matrix,
-        dict
+        boundaryCoeffs,
+        interfaces,
+        solverDict
     ),
     preconPtr_
     (
@@ -85,7 +89,13 @@ Foam::BlockBiCGStabSolver<Type>::solve
     Field<Type> p(x.size());
 
     // Calculate initial residual
-    matrix.Amul(p, x);
+    matrix.Amul
+    (
+        p,
+        x,
+        BlockLduSolver<Type>::boundaryCoeffs_,
+        BlockLduSolver<Type>::interfaces_
+    );
     Field<Type> r(b - p);
 
     solverPerf.initialResidual() = gSum(cmptMag(r))/norm;
@@ -138,7 +148,13 @@ Foam::BlockBiCGStabSolver<Type>::solve
             }
 
             preconPtr_->precondition(ph, p);
-            matrix.Amul(v, ph);
+            matrix.Amul
+            (
+                v,
+                ph,
+                BlockLduSolver<Type>::boundaryCoeffs_,
+                BlockLduSolver<Type>::interfaces_
+            );
             alpha = rho/gSumProd(rw, v);
 
             forAll (s, i)
@@ -147,7 +163,13 @@ Foam::BlockBiCGStabSolver<Type>::solve
             }
 
             preconPtr_->preconditionT(sh, s);
-            matrix.Amul(t, sh);
+            matrix.Amul
+            (
+                t,
+                sh,
+                BlockLduSolver<Type>::boundaryCoeffs_,
+                BlockLduSolver<Type>::interfaces_
+            );
             omega = gSumProd(t, s)/gSumProd(t, t);
 
             forAll (x, i)

src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockBiCGStab/BlockBiCGStabSolver.H

@@ -85,7 +85,9 @@ public:
         (
             const word& fieldName,
             const BlockLduMatrix<Type>& matrix,
-            const dictionary& dict
+            const FieldField<CoeffField, Type>& boundaryCoeffs,
+            const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
+            const dictionary& solverDict
         );
 
 

src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockCG/BlockCGSolver.C

@@ -37,14 +37,18 @@ Foam::BlockCGSolver<Type>::BlockCGSolver
 (
     const word& fieldName,
     const BlockLduMatrix<Type>& matrix,
-    const dictionary& dict
+    const FieldField<CoeffField, Type>& boundaryCoeffs,
+    const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
+    const dictionary& solverDict
 )
 :
     BlockIterativeSolver<Type>
     (
         fieldName,
         matrix,
-        dict
+        boundaryCoeffs,
+        interfaces,
+        solverDict
     ),
     preconPtr_
     (
@@ -84,7 +88,13 @@ typename Foam::BlockSolverPerformance<Type> Foam::BlockCGSolver<Type>::solve
     Field<Type> wA(x.size());
 
     // Calculate initial residual
-    matrix.Amul(wA, x);
+    matrix.Amul
+    (
+        wA,
+        x,
+        BlockLduSolver<Type>::boundaryCoeffs_,
+        BlockLduSolver<Type>::interfaces_
+    );
     Field<Type> rA(b - wA);
 
     solverPerf.initialResidual() = gSum(cmptMag(rA))/norm;
@@ -119,7 +129,13 @@ typename Foam::BlockSolverPerformance<Type> Foam::BlockCGSolver<Type>::solve
             }
 
             // Update preconditioner residual
-            matrix.Amul(wA, pA);
+            matrix.Amul
+            (
+                wA,
+                pA,
+                BlockLduSolver<Type>::boundaryCoeffs_,
+                BlockLduSolver<Type>::interfaces_
+            );
 
             wApA = gSumProd(wA, pA);
 

src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockCG/BlockCGSolver.H

@@ -85,7 +85,9 @@ public:
         (
             const word& fieldName,
             const BlockLduMatrix<Type>& matrix,
-            const dictionary& dict
+            const FieldField<CoeffField, Type>& boundaryCoeffs,
+            const typename BlockLduInterfaceFieldPtrsList<Type>::Type& interfaces,
+            const dictionary& solverDict
         );