Merge remote-tracking branch 'origin/hotfix/ggiSeparation+blockCoupledGgi' into nextRelease

Conflicts:
	src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.H

.gitignore

@@ -101,4 +101,12 @@ src/lduSolvers/amg/amgPolicy/samgPolicy.H
 src/lduSolvers/amg/amgPolicy/aamgPolicy.C
 src/lduSolvers/amg/amgPolicy/aamgPolicy.H
 
+# The following files are blacklisted because of a DMCA complaint by ANSYS.
+src/lduSolvers/tools/PriorityArray.C
+src/lduSolvers/tools/PriorityArray.H
+src/lduSolvers/amg/amgPolicy/samgPolicy.C
+src/lduSolvers/amg/amgPolicy/samgPolicy.H
+src/lduSolvers/amg/amgPolicy/aamgPolicy.C
+src/lduSolvers/amg/amgPolicy/aamgPolicy.H
+
 # end-of-file

src/OpenFOAM/interpolations/GGIInterpolation/GGIInterpolation.H

@@ -95,7 +95,7 @@ public:
         ClassName("GGIInterpolation");
 
         //- Quick reject names
-        static const NamedEnum<quickReject, 3> quickRejectNames_;
+        static const NamedEnum<quickReject, 4> quickRejectNames_;
 
 
     // Constructors

src/OpenFOAM/interpolations/GGIInterpolation/GGIInterpolationName.C

@@ -41,12 +41,15 @@ defineTypeNameAndDebug(Foam::GGIInterpolationName, 0);
 
 template<>
 const char*
-Foam::NamedEnum<Foam::GGIInterpolationName::quickReject, 3>::names[] =
+Foam::NamedEnum<Foam::GGIInterpolationName::quickReject, 4>::names[] =
 {
-    "3DDistance",
+    "distance3D",
     "AABB",
+    "bbOctree",
     "nSquared"
 };
 
+const Foam::NamedEnum<Foam::GGIInterpolationName::quickReject, 4>
+    Foam::GGIInterpolationName::quickRejectNames_;
 
 // ************************************************************************* //

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

@@ -64,15 +64,44 @@ void Foam::BlockLduMatrix<Type>::AmulCore
     const unallocLabelList& u = lduAddr().upperAddr();
     const unallocLabelList& l = lduAddr().lowerAddr();
 
+    
     const TypeCoeffField& Diag = this->diag();
     const TypeCoeffField& Upper = this->upper();
 
-    // Diagonal multiplication, no indirection
-    multiply(Ax, Diag, x);
-
     // Create multiplication function object
     typename BlockCoeff<Type>::multiply mult;
     
+    // AmulCore must be additive to account for initialisation step
+    // in ldu interfaces.  HJ, 6/Nov/2007
+    // Fixed by IC 19/Oct/2011
+    if (Diag.activeType() == blockCoeffBase::SCALAR)
+    {
+        const scalarTypeField& activeDiag = Diag.asScalar();
+
+        for (register label cellI = 0; cellI < u.size(); cellI++)
+        {
+            Ax[cellI] += mult(activeDiag[cellI], x[cellI]);
+        }
+    }
+    else if (Diag.activeType() == blockCoeffBase::LINEAR)
+    {
+        const linearTypeField& activeDiag = Diag.asLinear();
+
+        for (register label cellI = 0; cellI < u.size(); cellI++)
+        {
+            Ax[cellI] += mult(activeDiag[cellI], x[cellI]);
+        }
+    }
+    else if (Diag.activeType() == blockCoeffBase::SQUARE)
+    {
+        const squareTypeField& activeDiag = Diag.asSquare();
+
+        for (register label cellI = 0; cellI < u.size(); cellI++)
+        {
+            Ax[cellI] += mult(activeDiag[cellI], x[cellI]);
+        }
+    }
+
     // Lower multiplication
 
     if (symmetric())
@@ -209,12 +238,40 @@ void Foam::BlockLduMatrix<Type>::TmulCore
     const TypeCoeffField& Diag = this->diag();
     const TypeCoeffField& Upper = this->upper();
 
-    // Diagonal multiplication, no indirection
-    multiply(Tx, Diag, x);
-
     // Create multiplication function object
     typename BlockCoeff<Type>::multiply mult;
 
+    // AmulCore must be additive to account for initialisation step
+    // in ldu interfaces.  HJ, 6/Nov/2007
+    // Fixed by IC 19/Oct/2011
+    if (Diag.activeType() == blockCoeffBase::SCALAR)
+    {
+        const scalarTypeField& activeDiag = Diag.asScalar();
+
+        for (register label cellI = 0; cellI < u.size(); cellI++)
+        {
+            Tx[cellI] += mult(activeDiag[cellI], x[cellI]);
+        }
+    }
+    else if (Diag.activeType() == blockCoeffBase::LINEAR)
+    {
+        const linearTypeField& activeDiag = Diag.asLinear();
+
+        for (register label cellI = 0; cellI < u.size(); cellI++)
+        {
+            Tx[cellI] += mult(activeDiag[cellI], x[cellI]);
+        }
+    }
+    else if (Diag.activeType() == blockCoeffBase::SQUARE)
+    {
+        const squareTypeField& activeDiag = Diag.asSquare();
+
+        for (register label cellI = 0; cellI < u.size(); cellI++)
+        {
+            Tx[cellI] += mult(activeDiag[cellI].T(), x[cellI]);
+        }
+    }
+    
     // Upper multiplication
 
     if (Upper.activeType() == blockCoeffBase::SCALAR)

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

@@ -248,23 +248,31 @@ void BlockLduMatrix<scalar>::AmulCore
     const scalarField& x
 ) const
 {
-    const scalarField& Diag = diag();
-
-    for (register label rowI = 0; rowI < x.size(); rowI++)
-    {
-        Ax[rowI] = Diag[rowI]*x[rowI];
-    }
-
     // Note: pointer looping
     const label* const __restrict__ U = lduAddr().upperAddr().begin();
     const label* const __restrict__ L = lduAddr().lowerAddr().begin();
     
+    const scalar* const __restrict__ X = x.begin();
+    scalar* __restrict__ AX = Ax.begin();
+    
+    if (thereIsDiag())
+    {
+        const scalar* const __restrict__ diagPtr = diag().begin();
+
+        register const label nCells = diag().size();
+        for (register label cell=0; cell<nCells; cell++)
+        {
+            // AmulCore must be additive to account for initialisation step
+            // in ldu interfaces.  HJ, 6/Nov/2007
+            AX[cell] += diagPtr[cell]*X[cell];
+        }
+    }
+    
     if (symmetric())
+    {
+        if (thereIsUpper())
         {
             const scalar* const __restrict__ Upper = upper().begin();
-        const scalar* const __restrict__ X = x.begin();
-
-        scalar* __restrict__ AX = Ax.begin();
 
             for (register label coeffI = 0; coeffI < upper().size(); coeffI++)
             {
@@ -272,13 +280,21 @@ void BlockLduMatrix<scalar>::AmulCore
                 AX[L[coeffI]] += Upper[coeffI]*X[U[coeffI]];
             }
         }
+        else
+        {
+            const scalar* const __restrict__ Lower = lower().begin();
+
+            for (register label coeffI = 0; coeffI < upper().size(); coeffI++)
+            {
+                AX[U[coeffI]] += Lower[coeffI]*X[L[coeffI]];
+                AX[L[coeffI]] += Lower[coeffI]*X[U[coeffI]];
+            }
+        }
+    }
     else if (asymmetric())
     {
         const scalar* const __restrict__ Lower = lower().begin();
         const scalar* const __restrict__  Upper = upper().begin();
-        const scalar* const __restrict__ X = x.begin();
-
-        scalar* __restrict__ AX = Ax.begin();
 
         for (register label coeffI = 0; coeffI < upper().size(); coeffI++)
         {
@@ -296,23 +312,31 @@ void BlockLduMatrix<scalar>::TmulCore
     const scalarField& x
 ) const
 {
-    const scalarField& Diag = diag();
-
-    for (register label rowI = 0; rowI < x.size(); rowI++)
-    {
-        Tx[rowI] = Diag[rowI]*x[rowI];
-    }
-
     // Note: pointer looping
     const label* const __restrict__ U = lduAddr().upperAddr().begin();
     const label* const __restrict__ L = lduAddr().lowerAddr().begin();
     
+    const scalar* const __restrict__ X = x.begin();
+    scalar* __restrict__ TX = Tx.begin();
+    
+    if (thereIsDiag())
+    {
+        const scalar* const __restrict__ diagPtr = diag().begin();
+
+        register const label nCells = diag().size();
+        for (register label cell=0; cell<nCells; cell++)
+        {
+            // AmulCore must be additive to account for initialisation step
+            // in ldu interfaces.  HJ, 6/Nov/2007
+            TX[cell] += diagPtr[cell]*X[cell];
+        }
+    }
+
     if (symmetric())
+    {
+        if (thereIsUpper())
         {
             const scalar* const __restrict__ Upper = upper().begin();
-        const scalar* const __restrict__ X = x.begin();
-
-        scalar* __restrict__ TX = Tx.begin();
 
             for (register label coeffI = 0; coeffI < upper().size(); coeffI++)
             {
@@ -320,13 +344,21 @@ void BlockLduMatrix<scalar>::TmulCore
                 TX[L[coeffI]] += Upper[coeffI]*X[U[coeffI]];
             }
         }
+        else
+        {
+            const scalar* const __restrict__ Lower = lower().begin();
+
+            for (register label coeffI = 0; coeffI < upper().size(); coeffI++)
+            {
+                TX[U[coeffI]] += Lower[coeffI]*X[L[coeffI]];
+                TX[L[coeffI]] += Lower[coeffI]*X[U[coeffI]];
+            }
+        }
+    }
     else if (asymmetric())
     {
         const scalar* const __restrict__ Lower = lower().begin();
         const scalar* const __restrict__ Upper = upper().begin();
-        const scalar* const __restrict__ X = x.begin();
-
-        scalar* __restrict__ TX = Tx.begin();
 
         for (register label coeffI = 0; coeffI < upper().size(); coeffI++)
         {
@@ -368,9 +400,47 @@ tmp<scalarField> BlockLduMatrix<scalar>::H(const scalarField& x) const
         const label* const __restrict__ U = lduAddr().upperAddr().begin();
         const label* const __restrict__ L = lduAddr().lowerAddr().begin();
 
+        const scalar* const __restrict__ X = x.begin();
+        
+        if (symmetric())
+        {
+            if (thereIsUpper())
+            {
+                const scalar* const __restrict__ Upper = upper().begin();
+            
+                scalar* __restrict__ R = result.begin();
+        
+                for (register label coeffI = 0; coeffI < upper().size(); coeffI++)
+                {
+                    R[U[coeffI]] -= Upper[coeffI]*X[U[coeffI]];
+                }
+
+                for (register label coeffI = 0; coeffI < upper().size(); coeffI++)
+                {
+                    R[L[coeffI]] -= Upper[coeffI]*X[L[coeffI]];
+                }
+            }
+            else if (thereIsLower())
+            {
+                const scalar* const __restrict__ Lower = lower().begin();
+                
+                scalar* __restrict__ R = result.begin();
+
+                for (register label coeffI = 0; coeffI < upper().size(); coeffI++)
+                {
+                    R[U[coeffI]] -= Lower[coeffI]*X[U[coeffI]];
+                }
+
+                for (register label coeffI = 0; coeffI < upper().size(); coeffI++)
+                {
+                    R[L[coeffI]] -= Lower[coeffI]*X[L[coeffI]];
+                }
+            }
+        }
+        else
+        {
             const scalar* const __restrict__ Lower = lower().begin();
             const scalar* const __restrict__ Upper = upper().begin();
-        const scalar* const __restrict__ X = x.begin();
             
             scalar* __restrict__ R = result.begin();
 
@@ -384,6 +454,7 @@ tmp<scalarField> BlockLduMatrix<scalar>::H(const scalarField& x) const
                 R[L[coeffI]] -= Lower[coeffI]*X[L[coeffI]];
             }
         }
+    }
     
     return tresult;
 }
@@ -396,6 +467,43 @@ tmp<scalarField> BlockLduMatrix<scalar>::faceH(const scalarField& x) const
     scalarField& result = tresult();
 
     if (thereIsUpper() || thereIsLower())
+    {
+        if (symmetric())
+        {
+            if (thereIsUpper())
+            {
+                // Note: pointer looping
+                const label* const __restrict__ U = lduAddr().upperAddr().begin();
+                const label* const __restrict__ L = lduAddr().lowerAddr().begin();
+
+                const scalar* const __restrict__ Upper = upper().begin();
+                const scalar* const __restrict__ X = x.begin();
+
+                scalar* __restrict__ R = result.begin();
+
+                for (register label coeffI = 0; coeffI < upper().size(); coeffI++)
+                {
+                    R[coeffI] = Upper[coeffI]*(X[U[coeffI]] - X[L[coeffI]]);
+                }
+            }
+            else if (thereIsLower())
+            {
+                // Note: pointer looping
+                const label* const __restrict__ U = lduAddr().upperAddr().begin();
+                const label* const __restrict__ L = lduAddr().lowerAddr().begin();
+
+                const scalar* const __restrict__ Lower = lower().begin();
+                const scalar* const __restrict__ X = x.begin();
+
+                scalar* __restrict__ R = result.begin();
+
+                for (register label coeffI = 0; coeffI < upper().size(); coeffI++)
+                {
+                    R[coeffI] = Lower[coeffI]*(X[U[coeffI]] - X[L[coeffI]]);
+                }
+            }
+        }
+        else
         {
             // Note: pointer looping
             const label* const __restrict__ U = lduAddr().upperAddr().begin();
@@ -412,6 +520,7 @@ tmp<scalarField> BlockLduMatrix<scalar>::faceH(const scalarField& x) const
                 R[coeffI] = Upper[coeffI]*X[U[coeffI]] - Lower[coeffI]*X[L[coeffI]];
             }
         }
+    }
 
     return tresult;
 }

src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterfaceFields/ggiLduInterfaceField/ggiLduInterfaceField.C

@@ -69,4 +69,5 @@ void Foam::ggiLduInterfaceField::transformCoupleField
 }
 
 
+
 // ************************************************************************* //

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclicGgi/cyclicGgiPolyPatch.C

@@ -257,7 +257,7 @@ const Foam::cyclicGgiPolyPatch& Foam::cyclicGgiPolyPatch::cyclicShadow() const
 }
 
 
-void Foam::cyclicGgiPolyPatch::calcTransforms()
+void Foam::cyclicGgiPolyPatch::calcTransforms() const
 {
     if (active() && debug)
     {

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclicGgi/cyclicGgiPolyPatch.H

@@ -75,7 +75,7 @@ class cyclicGgiPolyPatch
 
         //- Calculate cyclic transforms (rotation and translation)
         //  Virtual over-ride for base GGI patch.  HJ, 14/Jan/2009
-        virtual void calcTransforms();
+        virtual void calcTransforms() const;
 
         //- Check definition: angles and offsets
         void checkDefinition() const;

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C

@@ -202,8 +202,7 @@ void Foam::ggiPolyPatch::calcPatchToPatch() const
                 0,             // Non-overlapping face tolerances
                 0,             // HJ, 24/Oct/2008
                 true,          // Rescale weighting factors.  Bug fix, MB.
-//                 ggiInterpolation::AABB
+                reject_        // Quick rejection algorithm, default BB_OCTREE
-                ggiInterpolation::BB_OCTREE  // Octree search, MB.
             );
 
         // Abort immediately if uncovered faces are present and the option
@@ -469,6 +468,7 @@ Foam::ggiPolyPatch::ggiPolyPatch
     shadowName_("initializeMe"),
     zoneName_("initializeMe"),
     bridgeOverlap_(false),
+    reject_(ggiZoneInterpolation::BB_OCTREE),
     shadowIndex_(-1),
     zoneIndex_(-1),
     patchToPatchPtr_(NULL),
@@ -490,13 +490,15 @@ Foam::ggiPolyPatch::ggiPolyPatch
     const polyBoundaryMesh& bm,
     const word& shadowName,
     const word& zoneName,
-    const bool bridgeOverlap
+    const bool bridgeOverlap,
+    const ggiZoneInterpolation::quickReject reject
 )
 :
     coupledPolyPatch(name, size, start, index, bm),
     shadowName_(shadowName),
     zoneName_(zoneName),
     bridgeOverlap_(bridgeOverlap),
+    reject_(reject),
     shadowIndex_(-1),
     zoneIndex_(-1),
     patchToPatchPtr_(NULL),
@@ -521,6 +523,7 @@ Foam::ggiPolyPatch::ggiPolyPatch
     shadowName_(dict.lookup("shadowPatch")),
     zoneName_(dict.lookup("zone")),
     bridgeOverlap_(dict.lookup("bridgeOverlap")),
+    reject_(ggiZoneInterpolation::BB_OCTREE),
     shadowIndex_(-1),
     zoneIndex_(-1),
     patchToPatchPtr_(NULL),
@@ -530,7 +533,15 @@ Foam::ggiPolyPatch::ggiPolyPatch
     localParallelPtr_(NULL),
     receiveAddrPtr_(NULL),
     sendAddrPtr_(NULL)
-{}
+{
+    if (dict.found("quickReject"))
+    {
+        reject_ = ggiZoneInterpolation::quickRejectNames_.read
+        (
+            dict.lookup("quickReject")
+        );
+    }
+}
 
 
 Foam::ggiPolyPatch::ggiPolyPatch
@@ -543,6 +554,7 @@ Foam::ggiPolyPatch::ggiPolyPatch
     shadowName_(pp.shadowName_),
     zoneName_(pp.zoneName_),
     bridgeOverlap_(pp.bridgeOverlap_),
+    reject_(pp.reject_),
     shadowIndex_(-1),
     zoneIndex_(-1),
     patchToPatchPtr_(NULL),
@@ -569,6 +581,7 @@ Foam::ggiPolyPatch::ggiPolyPatch
     shadowName_(pp.shadowName_),
     zoneName_(pp.zoneName_),
     bridgeOverlap_(pp.bridgeOverlap_),
+    reject_(pp.reject_),
     shadowIndex_(-1),
     zoneIndex_(-1),
     patchToPatchPtr_(NULL),
@@ -742,6 +755,11 @@ void Foam::ggiPolyPatch::initAddressing()
         // Calculate transforms for correct GGI cut
         calcTransforms();
         
+        if (master())
+        {
+            shadow().calcTransforms();
+        }
+
         // Force zone addressing and remote zone addressing
         // (uses GGI interpolator)
         zoneAddressing();
@@ -847,7 +865,7 @@ void Foam::ggiPolyPatch::updateMesh()
 }
 
 
-void Foam::ggiPolyPatch::calcTransforms()
+void Foam::ggiPolyPatch::calcTransforms() const
 {
     // Simplest mixing plane: no transform or separation.  HJ, 24/Oct/2008
     forwardT_.setSize(0);

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.H

@@ -73,6 +73,10 @@ class ggiPolyPatch
         //- Use bridging to fix overlap error in interpolation
         Switch bridgeOverlap_;
 
+        //- Quick reject algorithm
+        ggiZoneInterpolation::quickReject reject_;
+
+
         // Demand-driven data
 
             //- Shadow patch index.  Delayed evaluation for construction
@@ -214,7 +218,9 @@ public:
             const polyBoundaryMesh& bm,
             const word& shadowName,
             const word& zoneName,
-            const bool bridgeOverlap
+            const bool bridgeOverlap,
+            const ggiZoneInterpolation::quickReject
+                reject = ggiZoneInterpolation::BB_OCTREE
         );
 
         //- Construct from dictionary
@@ -307,6 +313,11 @@ public:
             //- Return interpolation face zone
             const faceZone& zone() const;
 
+        //- Is this the slave side?
+        bool slave() const
+        {
+            return !master();
+        }
 
         // Interpolation data
 

src/finiteVolume/fields/fvPatchFields/constraint/cyclicGgi/cyclicGgiFvPatchField.C

@@ -289,6 +289,35 @@ void cyclicGgiFvPatchField<Type>::updateInterfaceMatrix
 {}
 
 
+template<class Type>
+void cyclicGgiFvPatchField<Type>::initInterfaceMatrixUpdate
+(
+    const Field<Type>& psiInternal,
+    Field<Type>& result,
+    const BlockLduMatrix<Type>& m,
+    const CoeffField<Type>& coeffs,
+    const Pstream::commsTypes commsType
+) const
+{
+    notImplemented("void cyclicGgiFvPatchField::initInterfaceMatrixUpdate"
+        " for block-coupled matrices of tensor types"
+    )
+}
+
+
+template<class Type>
+void cyclicGgiFvPatchField<Type>::updateInterfaceMatrix
+(
+    const Field<Type>&,
+    Field<Type>&,
+    const BlockLduMatrix<Type>&,
+    const CoeffField<Type>&,
+    const Pstream::commsTypes commsType
+) const
+{
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/finiteVolume/fields/fvPatchFields/constraint/cyclicGgi/cyclicGgiFvPatchField.H

@@ -168,6 +168,25 @@ public:
                 const bool switchToLhs
             ) const;
 
+            //- Initialise neighbour matrix update
+            virtual void initInterfaceMatrixUpdate
+            (
+                const Field<Type>& psiInternal,
+                Field<Type>& result,
+                const BlockLduMatrix<Type>& m,
+                const CoeffField<Type>& coeffs,
+                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;
             
         //- GGI coupled interface functions
 

src/finiteVolume/fields/fvPatchFields/constraint/cyclicGgi/cyclicGgiFvPatchFields.C

@@ -37,12 +37,102 @@ Contributor:
 #include "cyclicGgiFvPatchFields.H"
 #include "addToRunTimeSelectionTable.H"
 #include "volFields.H"
+#include "scalarCoeffField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<>
+void cyclicGgiFvPatchField<scalar>::initInterfaceMatrixUpdate
+(
+    const Field<scalar>& psiInternal,
+    Field<scalar>& result,
+    const BlockLduMatrix<scalar>& m,
+    const CoeffField<scalar>& coeffs,
+    const Pstream::commsTypes commsType
+) const
+{
+    // Communication is allowed either before or after processor
+    // patch comms.  HJ, 11/Jul/2011
+
+    // Get shadow face-cells and assemble shadow field
+    const unallocLabelList& sfc = cyclicGgiPatch_.shadow().faceCells();
+
+    scalarField sField(sfc.size());
+
+    forAll (sField, i)
+    {
+        sField[i] = psiInternal[sfc[i]];
+    }
+
+    // Note: scalar interpolate does not get a transform, so this is safe
+    // HJ, 12/Jan/2009
+    scalarField pnf = cyclicGgiPatch_.interpolate(sField);
+    
+    // Multiply the field by coefficients and add into the result
+    const unallocLabelList& fc = cyclicGgiPatch_.faceCells();
+
+    forAll(fc, elemI)
+    {
+        result[fc[elemI]] -= coeffs[elemI]*pnf[elemI];
+    };
+}
+
+
+template<>
+void cyclicGgiFvPatchField<vector>::initInterfaceMatrixUpdate
+(
+    const Field<vector>& psiInternal,
+    Field<vector>& result,
+    const BlockLduMatrix<vector>& m,
+    const CoeffField<vector>& coeffs,
+    const Pstream::commsTypes commsType
+) const
+{
+    // Communication is allowed either before or after processor
+    // patch comms.  HJ, 11/Jul/2011
+
+    // Get shadow face-cells and assemble shadow field
+    const unallocLabelList& sfc = cyclicGgiPatch_.shadow().faceCells();
+
+    Field<vector> sField(sfc.size());
+
+    forAll (sField, i)
+    {
+        sField[i] = psiInternal[sfc[i]];
+    }
+
+    // Transformation is handled in interpolation.  HJ, 7/Jan/2009
+    Field<vector> pnf = cyclicGgiPatch_.interpolate(sField);
+
+    if (coeffs.activeType() == blockCoeffBase::SCALAR)                          
+    {                                                                           
+        pnf = coeffs.asScalar() * pnf;
+    }
+    else if (coeffs.activeType() == blockCoeffBase::LINEAR)                     
+    {                                                                           
+        pnf = cmptMultiply(coeffs.asLinear(), pnf);                             
+    }                                                                           
+    else if (coeffs.activeType() == blockCoeffBase::SQUARE)
+    {
+        pnf = coeffs.asSquare() & pnf;
+    }
+
+    // Multiply the field by coefficients and add into the result
+    const unallocLabelList& fc = cyclicGgiPatch_.faceCells();
+    
+    // Multiply the field by coefficients and add into the result
+    forAll(fc, elemI)
+    {
+        result[fc[elemI]] -= pnf[elemI];
+    }
+}
+
+
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 makePatchFields(cyclicGgi);

src/finiteVolume/fields/fvPatchFields/constraint/cyclicGgi/cyclicGgiFvPatchFields.H

@@ -51,6 +51,30 @@ SourceFiles
 namespace Foam
 {
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<>
+void cyclicGgiFvPatchField<scalar>::initInterfaceMatrixUpdate
+(
+    const Field<scalar>& psiInternal,
+    Field<scalar>& result,
+    const BlockLduMatrix<scalar>& m,
+    const CoeffField<scalar>& coeffs,
+    const Pstream::commsTypes commsType
+) const;
+
+
+template<>
+void cyclicGgiFvPatchField<vector>::initInterfaceMatrixUpdate
+(
+    const Field<vector>& psiInternal,
+    Field<vector>& result,
+    const BlockLduMatrix<vector>& m,
+    const CoeffField<vector>& coeffs,
+    const Pstream::commsTypes commsType
+) const;
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 makePatchTypeFieldTypedefs(cyclicGgi)