Serial Block GGI SAMG interface

src/foam/matrices/blockLduMatrix/BlockAMG/BlockAMGInterfaceFields/GGIBlockAMGInterfaceField/GGIBlockAMGInterfaceField.C

@@ -72,7 +72,7 @@ void Foam::GGIBlockAMGInterfaceField<Type>::agglomerateBlockType
     const scalarField& restrictWeights = ggiInterface_.restrictWeights();
 
     // Restrict coefficients
-    forAll(restrictAddressing, ffi)
+    forAll (restrictAddressing, ffi)
     {
         zoneCoarseCoeffs[restrictAddressing[ffi]] +=
             restrictWeights[ffi]*zoneFineCoeffs[fineAddressing[ffi]];
@@ -237,14 +237,14 @@ void Foam::GGIBlockAMGInterfaceField<Type>::updateInterfaceMatrix
 
     if (switchToLhs)
     {
-        forAll(faceCells, elemI)
+        forAll (faceCells, elemI)
         {
             result[faceCells[elemI]] += pnf[elemI];
         }
     }
     else
     {
-        forAll(faceCells, elemI)
+        forAll (faceCells, elemI)
         {
             result[faceCells[elemI]] -= pnf[elemI];
         }

src/foam/matrices/blockLduMatrix/BlockAMG/BlockSAMGInterfaceFields/BlockSAMGInterfaceField/blockVectorNSAMGInterfaceFields.C

@@ -34,6 +34,7 @@ Author
 
 #include "BlockSAMGInterfaceField.H"
 #include "ProcessorBlockSAMGInterfaceField.H"
+#include "GGIBlockSAMGInterfaceField.H"
 #include "VectorNFieldTypes.H"
 #include "ExpandTensorNField.H"
 
@@ -46,12 +47,15 @@ namespace Foam
 
 #define makeTemplateTypeNameAndDebug(type, Type, args...)                     \
                                                                               \
-typedef BlockSAMGInterfaceField<type > block##Type##SAMGInterfaceField;         \
+typedef BlockSAMGInterfaceField<type > block##Type##SAMGInterfaceField;       \
-defineNamedTemplateTypeNameAndDebug(block##Type##SAMGInterfaceField, 0);       \
+defineNamedTemplateTypeNameAndDebug(block##Type##SAMGInterfaceField, 0);      \
 defineTemplateRunTimeSelectionTable(block##Type##SAMGInterfaceField, lduInterface); \
                                                                               \
 typedef ProcessorBlockSAMGInterfaceField<type > block##Type##ProcessorSAMGInterfaceField;  \
 makeBlockSAMGInterfaceField(block##Type##SAMGInterfaceField, block##Type##ProcessorSAMGInterfaceField); \
+                                                                              \
+typedef GGIBlockSAMGInterfaceField<type > block##Type##GGISAMGInterfaceField; \
+makeBlockSAMGInterfaceField(block##Type##SAMGInterfaceField, block##Type##GGISAMGInterfaceField); \
 
 forAllVectorNTypes(makeTemplateTypeNameAndDebug);
 

src/foam/matrices/blockLduMatrix/BlockAMG/BlockSAMGInterfaceFields/GGIBlockSAMGInterfaceField/GGIBlockSAMGInterfaceField.C

@@ -0,0 +1,255 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend is distributed in the hope that it will be useful, but
+    WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+    General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "GGIBlockSAMGInterfaceField.H"
+#include "ggiLduInterfaceField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "blockLduMatrices.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+template<class Type2>
+void Foam::GGIBlockSAMGInterfaceField<Type>::selectBlockType
+(
+    Field<Type2>& coarseCoeffs,
+    const Foam::Field<Type2>& fineCoeffs
+) const
+{
+    // Note: reconsider better parallel communication here.
+    // Currently expanding to full zone size
+    // HJ, 16/Mar/2016
+
+    // Get fine interface
+    const ggiLduInterface& fineGgiInterface = ggiInterface_.fineGgiInterface();
+
+    // Reassemble fine coefficients to full fine zone size
+    // No need to initialise to zero, as only local coefficients
+    // are used.  HJ, 9/Jun/2016
+    Field<Type2> zoneFineCoeffs(fineGgiInterface.zoneSize());
+
+    const labelList& fineZa = fineGgiInterface.zoneAddressing();
+
+    forAll (fineZa, i)
+    {
+        zoneFineCoeffs[fineZa[i]] = fineCoeffs[i];
+    }
+
+    // Reduce zone data is not required: all coefficients are local
+    // HJ, 9/Jun/2016
+
+    Field<Type2> zoneCoarseCoeffs
+    (
+        ggiInterface_.zoneSize(),
+        pTraits<Type2>::zero
+    );
+
+    // Get addressing from the fine interface
+    const labelField& fineAddressing = ggiInterface_.fineAddressing();
+    const labelField& restrictAddressing = ggiInterface_.restrictAddressing();
+    const scalarField& restrictWeights = ggiInterface_.restrictWeights();
+
+    // Restrict coefficients
+    forAll (restrictAddressing, ffi)
+    {
+        zoneCoarseCoeffs[restrictAddressing[ffi]] +=
+            restrictWeights[ffi]*zoneFineCoeffs[fineAddressing[ffi]];
+    }
+
+    // Filter zone coefficients to local field
+    const labelList& za = ggiInterface_.zoneAddressing();
+
+    forAll (za, i)
+    {
+        coarseCoeffs[i] = zoneCoarseCoeffs[za[i]];
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::GGIBlockSAMGInterfaceField<Type>::GGIBlockSAMGInterfaceField
+(
+    const SAMGInterface& SAMGCp,
+    const BlockLduInterfaceField<Type>& fineInterfaceField
+)
+:
+    BlockSAMGInterfaceField<Type>(SAMGCp, fineInterfaceField),
+    ggiInterface_(refCast<const ggiSAMGInterface>(SAMGCp)),
+    doTransform_(false),
+    fieldTransferBuffer_()
+{
+    // If the interface based on a patch this must be taken care specially of
+    if (isA<GGIBlockLduInterfaceField<Type> >(fineInterfaceField))
+    {
+        const GGIBlockLduInterfaceField<Type>& p =
+            refCast<const GGIBlockLduInterfaceField<Type> >
+            (
+                fineInterfaceField
+            );
+
+        doTransform_ = p.doTransform();
+    }
+    else if (isA<ggiLduInterfaceField>(fineInterfaceField))
+    {
+        const ggiLduInterfaceField& p =
+            refCast<const ggiLduInterfaceField >(fineInterfaceField);
+
+        doTransform_ = p.doTransform();
+    }
+    else
+    {
+        FatalErrorIn("GGIBlockSAMGInterfaceField<Type> Constructor")
+            << "fineInterface must be of ggi type and either" << endl
+            << "    GGIBlockLduInterfaceField<Type> or " << endl
+            << "    ggiFvPatchField<Type> " << endl
+            << abort(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::GGIBlockSAMGInterfaceField<Type>::~GGIBlockSAMGInterfaceField()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp<Foam::CoeffField<Type> >
+Foam::GGIBlockSAMGInterfaceField<Type>::selectBlockCoeffs
+(
+    const Foam::CoeffField<Type>& fineCoeffs
+) const
+{
+    tmp<CoeffField<Type> > tcoarseCoeffs(new CoeffField<Type>(size()));
+    CoeffField<Type>& coarseCoeffs = tcoarseCoeffs();
+
+    typedef CoeffField<Type> TypeCoeffField;
+
+    typedef typename TypeCoeffField::linearTypeField linearTypeField;
+    typedef typename TypeCoeffField::squareTypeField squareTypeField;
+
+    // Added weights to account for non-integral matching
+    if (fineCoeffs.activeType() == blockCoeffBase::SQUARE)
+    {
+        squareTypeField& activeCoarseCoeffs = coarseCoeffs.asSquare();
+        const squareTypeField& activeFineCoeffs = fineCoeffs.asSquare();
+
+        this->selectBlockType(activeCoarseCoeffs, activeFineCoeffs);
+    }
+    else if (fineCoeffs.activeType() == blockCoeffBase::LINEAR)
+    {
+        linearTypeField& activeCoarseCoeffs = coarseCoeffs.asLinear();
+        const linearTypeField& activeFineCoeffs = fineCoeffs.asLinear();
+
+        this->selectBlockType(activeCoarseCoeffs, activeFineCoeffs);
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "Foam::tmp<Foam::CoeffField<Type> >\n"
+            "Foam::GGIBlockSAMGInterfaceField<Type>::selectBlockCoeffs\n"
+            "(\n"
+            "    const Foam::CoeffField<Type>& fineCoeffs\n"
+            ") const"
+        )   << "Scalar type agglomeration currently not handled"
+            << abort(FatalError);
+    }
+
+    return tcoarseCoeffs;
+}
+
+
+template<class Type>
+void Foam::GGIBlockSAMGInterfaceField<Type>::initInterfaceMatrixUpdate
+(
+    const Field<Type>& psiInternal,
+    Field<Type>&,
+    const BlockLduMatrix<Type>&,
+    const CoeffField<Type>&,
+    const Pstream::commsTypes,
+    const bool switchToLhs
+) const
+{
+    // This must have a reduce in it.  HJ, 15/May/2009
+    Field<Type> pif = ggiInterface_.interfaceInternalField(psiInternal);
+
+    fieldTransferBuffer_ = ggiInterface_.fastReduce(pif);
+}
+
+
+template<class Type>
+void Foam::GGIBlockSAMGInterfaceField<Type>::updateInterfaceMatrix
+(
+    const Field<Type>& psiInternal,
+    Field<Type>& result,
+    const BlockLduMatrix<Type>& matrix,
+    const CoeffField<Type>& coeffs,
+    const Pstream::commsTypes commsType,
+    const bool switchToLhs
+) const
+{
+    // Get interface from shadow
+    const GGIBlockSAMGInterfaceField<Type>& shadowInterface =
+        refCast<const GGIBlockSAMGInterfaceField<Type> >
+        (
+            matrix.interfaces()[ggiInterface_.shadowIndex()]
+        );
+
+    Field<Type> pnf = shadowInterface.fieldTransferBuffer();
+
+    // Complex (VectorN) transformation happens here.
+    // HJ, 17/Feb/2016
+//     transformCoupleField(pnf, cmpt);
+
+    // Multiply neighbour field with coeffs and re-use pnf for result
+    // of multiplication
+    multiply(pnf, coeffs, pnf);
+
+    const unallocLabelList& faceCells = ggiInterface_.faceCells();
+
+    if (switchToLhs)
+    {
+        forAll (faceCells, elemI)
+        {
+            result[faceCells[elemI]] += pnf[elemI];
+        }
+    }
+    else
+    {
+        forAll (faceCells, elemI)
+        {
+            result[faceCells[elemI]] -= pnf[elemI];
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/foam/matrices/blockLduMatrix/BlockAMG/BlockSAMGInterfaceFields/GGIBlockSAMGInterfaceField/GGIBlockSAMGInterfaceField.H

@@ -0,0 +1,224 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend is distributed in the hope that it will be useful, but
+    WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+    General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::GGIBlockSAMGInterfaceField
+
+Description
+    AMG selected GGI interface field.
+
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved
+
+SourceFiles
+    GGIBlockSAMGInterfaceField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef GGIBlockSAMGInterfaceField_H
+#define GGIBlockSAMGInterfaceField_H
+
+#include "BlockSAMGInterfaceField.H"
+#include "ggiSAMGInterface.H"
+#include "GGIBlockLduInterfaceField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class GGIBlockSAMGInterfaceField Declaration
+\*---------------------------------------------------------------------------*/
+
+template <class Type>
+class GGIBlockSAMGInterfaceField
+:
+    public BlockSAMGInterfaceField<Type>,
+    public GGIBlockLduInterfaceField<Type>
+{
+    // Private data
+
+        //- Local reference cast into the ggi interface
+        const ggiSAMGInterface& ggiInterface_;
+
+        //- Is the transform required
+        bool doTransform_;
+
+        //- Field transfer buffer
+        mutable Field<Type> fieldTransferBuffer_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        GGIBlockSAMGInterfaceField(const GGIBlockSAMGInterfaceField&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const GGIBlockSAMGInterfaceField&);
+
+
+        //- Select type
+        template<class Type2>
+        void selectBlockType
+        (
+            Field<Type2>& coarseCoeffs,
+            const Foam::Field<Type2>& fineCoeffs
+        ) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("ggi");
+
+
+    // Constructors
+
+        //- Construct from SAMG interface and fine level interface field
+        GGIBlockSAMGInterfaceField
+        (
+            const SAMGInterface& SAMGCp,
+            const BlockLduInterfaceField<Type>& fineInterfaceField
+        );
+
+
+    //- Destructor
+    virtual ~GGIBlockSAMGInterfaceField();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return size
+            label size() const
+            {
+                return ggiInterface_.size();
+            }
+
+
+        // Agglomeration
+
+            // Klas Jareteg: 2013-02-06. Moved by HJ, 16/Mar/2016
+            //- Agglomerating for the CoeffField fine-level coefficients
+            virtual tmp<CoeffField<Type> > selectBlockCoeffs
+            (
+                const CoeffField<Type>& fineCoeffs
+            ) const;
+
+
+        // Block coupled interface matrix update
+
+            //- Transform given patch component field
+            virtual void transformCoupleField
+            (
+                scalarField& f,
+                const direction cmpt
+            ) const
+            {
+                GGIBlockLduInterfaceField<Type>::transformCoupleField
+                (
+                    f,
+                    cmpt
+                );
+            }
+
+            //- Transform neighbour field
+            virtual void transformCoupleField
+            (
+                Field<Type>& f
+            ) const
+            {
+                GGIBlockLduInterfaceField<Type>::transformCoupleField(f);
+            }
+
+            //- Initialise neighbour matrix update
+            virtual void initInterfaceMatrixUpdate
+            (
+                const Field<Type>&,
+                Field<Type>&,
+                const BlockLduMatrix<Type>&,
+                const CoeffField<Type>&,
+                const Pstream::commsTypes commsType,
+                const bool switchToLhs
+            ) 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 bool switchToLhs
+            ) const;
+
+
+        //- Ggi interface functions
+
+            //- Does the interface field perform the transfromation
+            virtual bool doTransform() const
+            {
+                return doTransform_;
+            }
+
+            //- Return face transformation tensor
+            virtual const tensorField& forwardT() const
+            {
+                return ggiInterface_.forwardT();
+            }
+
+            //- Return neighbour-cell transformation tensor
+            virtual const tensorField& reverseT() const
+            {
+                return ggiInterface_.reverseT();
+            }
+
+
+        // Transfer buffer access
+
+            //- Return contents of the field transfer buffer
+            const Field<Type>& fieldTransferBuffer() const
+            {
+                return fieldTransferBuffer_;
+            }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "GGIBlockSAMGInterfaceField.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/foam/matrices/blockLduMatrix/BlockAMG/BlockSAMGInterfaceFields/GGIBlockSAMGInterfaceField/GGIBlockSAMGInterfaceFields.C

@@ -0,0 +1,42 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend is distributed in the hope that it will be useful, but
+    WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+    General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "GGIBlockSAMGInterfaceFields.H"
+#include "blockSAMGInterfaceFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+makeBlockSAMGInterfaceFields(GGIBlockSAMGInterfaceField);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/foam/matrices/blockLduMatrix/BlockAMG/BlockSAMGInterfaceFields/GGIBlockSAMGInterfaceField/GGIBlockSAMGInterfaceFields.H

@@ -0,0 +1,63 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend is distributed in the hope that it will be useful, but
+    WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+    General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    GGIBlockSAMGInterfaceField
+
+Description
+    Typedefs for block coefficient ggi SAMG interface fields
+
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved.
+
+SourceFiles
+    GGIBlockSAMGInterfaceFields.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef GGIBlockSAMGInterfaceFields_H
+#define GGIBlockSAMGInterfaceFields_H
+
+#include "GGIBlockSAMGInterfaceField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+typedef GGIBlockSAMGInterfaceField<scalar> GGIBlockSAMGInterfaceFieldScalar;
+typedef GGIBlockSAMGInterfaceField<vector> GGIBlockSAMGInterfaceFieldVector;
+typedef GGIBlockSAMGInterfaceField<tensor> GGIBlockSAMGInterfaceFieldTensor;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/foam/matrices/lduMatrix/solvers/AMG/interfaces/SAMGInterfaces/ggiSAMGInterface/ggiSAMGInterface.C

@@ -282,6 +282,20 @@ Foam::ggiSAMGInterface::ggiSAMGInterface
     // Note: local addressing contains only local faces
     const labelList& fineZa =  fineGgiInterface_.zoneAddressing();
 
+    // Expand master prolongation to zone
+    crMatrix masterExpandProlongation(interfaceProlongation);
+
+    //HJ, HERE: expand master prolongation to zone
+    // Note: master is now the size of local zone
+    if (!fineGgiInterface_.localParallel())
+    {
+        // Not line this: expand without communication HJ, HERE
+        // fineGgiInterface_.expandCrMatrixToZone
+        // (
+        //     masterExpandProlongation
+        // );
+    }
+
     // Create crMatrix for neighbour faces.  Note: expandCrMatrixToZone will
     // expand the matrix to zone size, including communications.
     // Faces which are not used locally will be marked by empty rows
@@ -298,10 +312,6 @@ Foam::ggiSAMGInterface::ggiSAMGInterface
         );
     }
 
-    // 
-    const crAddressing& prolongationCr = interfaceProlongation.crAddr();
-    const crAddressing& nbrExpandCr = nbrExpandProlongation.crAddr();
-
     // Create addressing for neighbour processors.  Note: expandAddrToZone will
     // expand the addressing to zone size.  HJ, 13/Jun/2016
     labelField neighbourExpandProc
@@ -369,6 +379,23 @@ Foam::ggiSAMGInterface::ggiSAMGInterface
     // Count the number of agglomeration pairs
     label nAgglomPairs = 0;
 
+    // Switching prolongation matrices
+    const crMatrix* masterP = NULL;
+    const crMatrix* neighbourP = NULL;
+
+    if (fineGgiInterface_.master())
+    {
+        // Grab prolongation matrix
+        masterP = &masterExpandProlongation;
+        neighbourP = &nbrExpandProlongation;
+    }
+    else
+    {
+        // Grab prolongation matrix
+        masterP = &nbrExpandProlongation;
+        neighbourP = &masterExpandProlongation;
+    }
+
     // On the fine level, addressing is made in a labelListList
     if (fineGgiInterface_.fineLevel())
     {
@@ -455,32 +482,31 @@ Foam::ggiSAMGInterface::ggiSAMGInterface
                 // restriction - with included weights from GGI
                 for
                 (
-                    label indexR = prolongationCr.rowStart()[curSide];
+                    label indexR = masterP->crAddr().rowStart()[curSide];
-                    indexR < prolongationCr.rowStart()[curSide + 1];
+                    indexR < masterP->crAddr().rowStart()[curSide + 1];
                     indexR++
                 )
                 {
                     // Grab weight from restriction
-                    scalar rWeight = interfaceProlongation.coeffs()[indexR];
+                    scalar rWeight = masterP->coeffs()[indexR];
 
-                    // HJ, replace nbrInterfaceProlongation with nbrExpandProlongation
                     for
                     (
-                        label indexP = nbrExpandCr.rowStart()[nbrSide];
+                        label indexP = neighbourP->crAddr().rowStart()[nbrSide];
-                        indexP < nbrExpandCr.rowStart()[nbrSide + 1];
+                        indexP < neighbourP->crAddr().rowStart()[nbrSide + 1];
                         indexP++
                     )
                     {
                         // Grab weight from prolongation
-                        scalar pWeight = nbrInterfaceProlongation.coeffs()[indexP];
+                        scalar pWeight = neighbourP->coeffs()[indexP];
 
                         // Code in the current master and slave - used for
                         // identifying the face
                         curMaster =
-                            prolongationCr.column()[indexR]
+                            masterP->crAddr().column()[indexR]
                           + procOffset*curMasterProc;
                         curSlave =
-                            nbrExpandCr.column()[indexP]
+                            neighbourP->crAddr().column()[indexP]
                           + procOffset*curSlaveProc;
 
                         if (neighboursTable.found(curMaster))
@@ -520,7 +546,10 @@ Foam::ggiSAMGInterface::ggiSAMGInterface
                                     nbrFound = true;
                                     curFaceFaces[curNbrI].append(ffI);
                                     curFaceFaceNbrs[curNbrI].append(nbrI);
-                                    curFaceWeights[curNbrI].append(curNW*pWeight*rWeight);
+                                    curFaceWeights[curNbrI].append
+                                    (
+                                        curNW*pWeight*rWeight
+                                    );
 
                                     // New agglomeration pair found in already
                                     // existing pair
@@ -609,16 +638,6 @@ Foam::ggiSAMGInterface::ggiSAMGInterface
                 }
             } // end for all current neighbours
         } // end for all fine faces
-        if (fineGgiInterface_.master())
-        {
-            Info<< "MASTER: ";
-        }
-        else
-        {
-            Info<< "SLAVE: ";
-        }
-        Info<< "Done fine level, 1. Created " << nAgglomPairs << " pairs and "
-            << nCoarseFaces << " faces" << endl;
     }
 //------------------------------------------------------------------------------
 //                          FINE LEVEL - no GGI weights!
@@ -685,30 +704,29 @@ Foam::ggiSAMGInterface::ggiSAMGInterface
             // restriction - with included weights from GGI
             for
             (
-                label indexR = prolongationCr.rowStart()[curSide];
+                label indexR = masterP->crAddr().rowStart()[curSide];
-                indexR < prolongationCr.rowStart()[curSide + 1];
+                indexR < masterP->crAddr().rowStart()[curSide + 1];
                 indexR++
             )
             {
                 // Grab weight from restriction
-                scalar rWeight = interfaceProlongation.coeffs()[indexR];
+                scalar rWeight = masterP->coeffs()[indexR];
 
-                // HJ, replace nbrInterfaceProlongation with nbrExpandProlongation
                 for
                 (
-                    label indexP = nbrExpandCr.rowStart()[nbrSide];
+                    label indexP = neighbourP->crAddr().rowStart()[nbrSide];
-                    indexP < nbrExpandCr.rowStart()[nbrSide + 1];
+                    indexP < neighbourP->crAddr().rowStart()[nbrSide + 1];
                     indexP++
                 )
                 {
                     // Grab weight from prolongation
-                    scalar pWeight = nbrInterfaceProlongation.coeffs()[indexP];
+                    scalar pWeight = neighbourP->coeffs()[indexP];
 
                     // Code in the current master and slave - used for
                     // identifying the face
-                    curMaster = prolongationCr.column()[indexR]
+                    curMaster = masterP->crAddr().column()[indexR]
                         + procOffset*curMasterProc;
-                    curSlave = nbrExpandCr.column()[indexP]
+                    curSlave = neighbourP->crAddr().column()[indexP]
                         + procOffset*curSlaveProc;
 
                     if (neighboursTable.found(curMaster))
@@ -835,8 +853,6 @@ Foam::ggiSAMGInterface::ggiSAMGInterface
                 }
             }
         } // end for all fine faces
-        Info<< "Done coarse level, 1. Created " << nAgglomPairs << " pairs and "
-            << nCoarseFaces << " faces" << endl;
     } // end of else in fine level (coarse level)
 
     // Since only local faces are analysed, lists can now be resized
@@ -894,7 +910,7 @@ Foam::ggiSAMGInterface::ggiSAMGInterface
     // Sort makes sure the order is identical on both sides.
     // HJ, 20/Feb/2009 and 6/Jun/2016
     sort(contents);
-    Info<< "START MATRIX ASSEMBLY" << endl;
+
     // Note: Restriction is done on master side only because this is where
     // the local zone is created.  HJ, 1/Aug/2016
     if (master())
@@ -983,8 +999,6 @@ Foam::ggiSAMGInterface::ggiSAMGInterface
                 nProcFaces++;
             }
         }
-        Info<< "MASTER ASSEMBLY: Created " << nAgglomPairs << " pairs and "
-            << nProcFaces << " master faces" << endl;
 
         // No need to resize arrays only local faces are used
         // HJ, 1/Aug/2016
@@ -1137,11 +1151,7 @@ Foam::ggiSAMGInterface::ggiSAMGInterface
                 nProcFaces++;
             }
         }
-        Info<< "SLAVE ASSEMBLY: Created " << nAgglomPairs << " pairs and "
-            << nProcFaces << " slave faces" << endl;
-
     }
-    Info<< "END MATRIX ASSEMBLY" << endl;
 }
 
 
@@ -1337,7 +1347,6 @@ void Foam::ggiSAMGInterface::expandAddrToZone(labelField& lf) const
 
 void Foam::ggiSAMGInterface::expandCrMatrixToZone(crMatrix&) const
 {
-    notImplemented("expandCrMatrixToZone");
     // Code missing: collapse crMatrices into a zone crMatrix
     if (!localParallel())
     {
@@ -1353,9 +1362,8 @@ void Foam::ggiSAMGInterface::initProlongationTransfer
     const crMatrix& filteredP
 ) const
 {
-    // Send prolongation matrix, using IOstream operators
+    // crMatrix transfer is local without global reduction
-    //OPstream toNbr(Pstream::blocking, neighbProcNo());
+    crMatrixTransferBuffer_ = filteredP;
-    //toNbr<< filteredP;
 }
 
 
@@ -1366,11 +1374,13 @@ Foam::ggiSAMGInterface::prolongationTransfer
     const crMatrix& filteredP
 ) const
 {
-    //IPstream fromNbr(Pstream::blocking, neighbProcNo());
+    autoPtr<crMatrix> tnbrP
-
+    (
-    autoPtr<crMatrix> tnbrP(new crMatrix(5,5, labelList(5,1)));
+        new crMatrix(shadowInterface().crMatrixTransferBuffer())
+    );
 
     return tnbrP;
 }
 
+
 // ************************************************************************* //