Implicit treatment of coupled boundary conditions

src/OpenFOAM/Make/files

@@ -431,6 +431,7 @@ $(lduAddressing)/lduAddressing.C
 $(lduAddressing)/lduInterface/lduInterface.C
 $(lduAddressing)/lduInterface/processorLduInterface.C
 $(lduAddressing)/lduInterface/cyclicLduInterface.C
+$(lduAddressing)/lduInterface/lduPrimitiveProcessorInterface.C
 
 lduInterfaceFields = $(lduAddressing)/lduInterfaceFields
 $(lduInterfaceFields)/lduInterfaceField/lduInterfaceField.C

src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduInterfaceField/LduInterfaceField.H

@@ -48,8 +48,7 @@ SourceFiles
 namespace Foam
 {
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
+// Forward declarations
 class lduMatrix;
 
 /*---------------------------------------------------------------------------*\
@@ -101,6 +100,8 @@ public:
             (
                 Field<Type>&,
                 const bool add,
+                const lduAddressing&,
+                const label interfacei,
                 const Field<Type>&,
                 const scalarField&,
                 const Pstream::commsTypes commsType
@@ -115,6 +116,8 @@ public:
             (
                 Field<Type>&,
                 const bool add,
+                const lduAddressing&,
+                const label interfacei,
                 const Field<Type>&,
                 const scalarField&,
                 const Pstream::commsTypes commsType

src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixUpdateMatrixInterfaces.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -53,6 +54,8 @@ void Foam::LduMatrix<Type, DType, LUType>::initMatrixInterfaces
                 (
                     result,
                     add,
+                    lduMesh_.lduAddr(),
+                    interfacei,
                     psiif,
                     interfaceCoeffs[interfacei],
                     //Amultiplier<Type, LUType>(interfaceCoeffs[interfacei]),
@@ -80,6 +83,8 @@ void Foam::LduMatrix<Type, DType, LUType>::initMatrixInterfaces
                 (
                     result,
                     add,
+                    lduMesh_.lduAddr(),
+                    interfacei,
                     psiif,
                     interfaceCoeffs[interfacei],
                     //Amultiplier<Type, LUType>(interfaceCoeffs[interfacei]),
@@ -128,6 +133,8 @@ void Foam::LduMatrix<Type, DType, LUType>::updateMatrixInterfaces
                 (
                     result,
                     add,
+                    lduMesh_.lduAddr(),
+                    interfacei,
                     psiif,
                     interfaceCoeffs[interfacei],
                     //Amultiplier<Type, LUType>(interfaceCoeffs[interfacei]),
@@ -153,6 +160,8 @@ void Foam::LduMatrix<Type, DType, LUType>::updateMatrixInterfaces
                     (
                         result,
                         add,
+                        lduMesh_.lduAddr(),
+                        interfacei,
                         psiif,
                         interfaceCoeffs[interfacei],
                       //Amultiplier<Type, LUType>(interfaceCoeffs[interfacei]),
@@ -165,6 +174,8 @@ void Foam::LduMatrix<Type, DType, LUType>::updateMatrixInterfaces
                     (
                         result,
                         add,
+                        lduMesh_.lduAddr(),
+                        interfacei,
                         psiif,
                         interfaceCoeffs[interfacei],
                       //Amultiplier<Type, LUType>(interfaceCoeffs[interfacei]),
@@ -189,6 +200,8 @@ void Foam::LduMatrix<Type, DType, LUType>::updateMatrixInterfaces
                 (
                     result,
                     add,
+                    lduMesh_.lduAddr(),
+                    interfacei,
                     psiif,
                     interfaceCoeffs[interfacei],
                     //Amultiplier<Type, LUType>(interfaceCoeffs[interfacei]),

src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduInterface.H

@@ -72,6 +72,7 @@ public:
 
     // Member Functions
 
+
         // Access
 
             //- Return faceCell addressing
@@ -87,6 +88,14 @@ public:
                 const labelUList& internalData
             ) const = 0;
 
+            //- Return the values of the given internal data adjacent to
+            //- the interface as a field using faceCell mapping
+            virtual tmp<labelField> interfaceInternalField
+            (
+                const labelUList& internalData,
+                const labelUList& faceCells
+            ) const = 0;
+
             //- Initialise transfer of internal field adjacent to the interface
             virtual void initInternalFieldTransfer
             (
@@ -95,12 +104,23 @@ public:
             ) const
             {}
 
+            //- Initialise transfer of internal field adjacent to the interface
+            //- using faceCells mapping
+            virtual void initInternalFieldTransfer
+            (
+                const Pstream::commsTypes commsType,
+                const labelUList& iF,
+                const labelUList& faceCells
+            ) const
+            {}
+
             //- Transfer and return internal field adjacent to the interface
             virtual tmp<labelField> internalFieldTransfer
             (
                 const Pstream::commsTypes commsType,
                 const labelUList& iF
             ) const = 0;
+
 };
 
 

src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduPrimitiveProcessorInterface.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2017 OpenCFD Ltd.
+    Copyright (C) 2016-2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -65,11 +65,23 @@ Foam::lduPrimitiveProcessorInterface::interfaceInternalField
     const labelUList& internalData
 ) const
 {
-    tmp<labelField> tfld(new labelField(faceCells_.size()));
-    labelField& fld = tfld.ref();
-    forAll(faceCells_, i)
+    return interfaceInternalField(internalData, faceCells_);
+}
+
+
+Foam::tmp<Foam::labelField>
+Foam::lduPrimitiveProcessorInterface::interfaceInternalField
+(
+    const labelUList& internalData,
+    const labelUList& faceCells
+) const
+{
+    auto tfld = tmp<labelField>::New(faceCells.size());
+    auto& fld = tfld.ref();
+
+    forAll(faceCells, i)
     {
-        fld[i] = internalData[faceCells_[i]];
+        fld[i] = internalData[faceCells[i]];
     }
     return tfld;
 }
@@ -85,6 +97,21 @@ void Foam::lduPrimitiveProcessorInterface::initInternalFieldTransfer
 }
 
 
+void Foam::lduPrimitiveProcessorInterface::initInternalFieldTransfer
+(
+    const Pstream::commsTypes commsType,
+    const labelUList& iF,
+    const labelUList& faceCells
+) const
+{
+    processorLduInterface::send
+    (
+        commsType,
+        interfaceInternalField(iF, faceCells)()
+    );
+}
+
+
 Foam::tmp<Foam::labelField>
 Foam::lduPrimitiveProcessorInterface::internalFieldTransfer
 (
@@ -95,5 +122,4 @@ Foam::lduPrimitiveProcessorInterface::internalFieldTransfer
     return processorLduInterface::receive<label>(commsType, faceCells_.size());
 }
 
-
 // ************************************************************************* //

src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduPrimitiveProcessorInterface.H

@@ -56,7 +56,7 @@ class lduPrimitiveProcessorInterface
     public lduInterface,
     public processorLduInterface
 {
-    // Private data
+    // Private Data
 
         const labelList faceCells_;
 
@@ -126,12 +126,20 @@ public:
         // Interface transfer functions
 
             //- Return the values of the given internal data adjacent to
-            //  the interface as a field
+            //- the interface as a field
             virtual tmp<labelField> interfaceInternalField
             (
                 const labelUList& internalData
             ) const;
 
+            //- Return the values of the given internal data adjacent to
+            //- the interface as a field using faceCells mapping
+            virtual tmp<labelField> interfaceInternalField
+            (
+                const labelUList& internalData,
+                const labelUList& faceCells
+            ) const;
+
             //- Initialise neighbour field transfer
             virtual void initInternalFieldTransfer
             (
@@ -139,6 +147,14 @@ public:
                 const labelUList& iF
             ) const;
 
+            //- Initialise neighbour field transfer using faceCells
+            virtual void initInternalFieldTransfer
+            (
+                const Pstream::commsTypes commsType,
+                const labelUList& iF,
+                const labelUList& faceCells
+            ) const;
+
             //- Transfer and return internal field adjacent to the interface
             virtual tmp<labelField> internalFieldTransfer
             (

src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterfaceFields/lduInterfaceField/lduInterfaceField.H

@@ -42,14 +42,14 @@ SourceFiles
 #include "lduInterface.H"
 #include "primitiveFieldsFwd.H"
 #include "Pstream.H"
+#include "lduAddressing.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
+// Forward declarations
 class lduMatrix;
 
 /*---------------------------------------------------------------------------*\
@@ -58,7 +58,7 @@ class lduMatrix;
 
 class lduInterfaceField
 {
-    // Private data
+    // Private Data
 
         //- Reference to the coupled patch this field is defined for
         const lduInterface& interface_;
@@ -140,6 +140,8 @@ public:
             (
                 solveScalarField& result,
                 const bool add,
+                const lduAddressing&,
+                const label interfacei,
                 const solveScalarField& psiInternal,
                 const scalarField& coeffs,
                 const direction cmpt,
@@ -153,6 +155,8 @@ public:
             (
                 solveScalarField& result,
                 const bool add,
+                const lduAddressing&,
+                const label interfacei,
                 const solveScalarField& psiInternal,
                 const scalarField& coeffs,
                 const direction cmpt,
@@ -165,6 +169,7 @@ public:
             (
                 Field<Type>& result,
                 const bool add,
+                const labelUList& faceCells,
                 const scalarField& coeffs,
                 const Field<Type>& vals
             ) const;

src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterfaceFields/lduInterfaceField/lduInterfaceFieldTemplates.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017 OpenCFD Ltd.
+    Copyright (C) 2017-2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,24 +32,23 @@ void Foam::lduInterfaceField::addToInternalField
 (
     Field<Type>& result,
     const bool add,
+    const labelUList& faceCells,
     const scalarField& coeffs,
     const Field<Type>& vals
 ) const
 {
-    const labelUList& faceCells = this->interface().faceCells();
-
     if (add)
     {
-        forAll(faceCells, elemI)
+        forAll(faceCells, elemi)
         {
-            result[faceCells[elemI]] += coeffs[elemI]*vals[elemI];
+            result[faceCells[elemi]] += coeffs[elemi]*vals[elemi];
         }
     }
     else
     {
-        forAll(faceCells, elemI)
+        forAll(faceCells, elemi)
         {
-            result[faceCells[elemI]] -= coeffs[elemI]*vals[elemI];
+            result[faceCells[elemi]] -= coeffs[elemi]*vals[elemi];
         }
     }
 }

src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C

@@ -325,7 +325,7 @@ void Foam::lduMatrix::setResidualField
     const bool initial
 ) const
 {
-    if (!lduMesh_.hasDb())
+    if (!mesh().hasDb())
     {
         return;
     }
@@ -341,12 +341,12 @@ void Foam::lduMatrix::setResidualField
     }
 
     scalarIOField* residualPtr =
-        lduMesh_.thisDb().getObjectPtr<scalarIOField>(lookupName);
+        mesh().thisDb().getObjectPtr<scalarIOField>(lookupName);
 
     if (residualPtr)
     {
         const IOdictionary* dataPtr =
-            lduMesh_.thisDb().findObject<IOdictionary>("data");
+            mesh().thisDb().findObject<IOdictionary>("data");
 
         if (dataPtr)
         {

src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H

@@ -86,7 +86,8 @@ class lduMatrix
     // private data
 
         //- LDU mesh reference
-        const lduMesh& lduMesh_;
+        //const lduMesh& lduMesh_;
+        std::reference_wrapper<const lduMesh> lduMesh_;
 
         //- Coefficients (not including interfaces)
         scalarField *lowerPtr_, *diagPtr_, *upperPtr_;
@@ -569,10 +570,16 @@ public:
                 return lduMesh_;
             }
 
+            //- Set the LDU mesh containing the addressing is obtained
+            void setLduMesh(const lduMesh& m)
+            {
+                lduMesh_ = m;
+            }
+
             //- Return the LDU addressing
             const lduAddressing& lduAddr() const
             {
-                return lduMesh_.lduAddr();
+                return mesh().lduAddr();
             }
 
             //- Return the patch evaluation schedule

src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixSolver.C

@@ -204,13 +204,13 @@ Foam::solveScalarField::cmptType Foam::lduMatrix::solver::normFactor
     // --- Calculate A dot reference value of psi
     matrix_.sumA(tmpField, interfaceBouCoeffs_, interfaces_);
 
-    tmpField *= gAverage(psi, matrix_.lduMesh_.comm());
+    tmpField *= gAverage(psi, matrix_.mesh().comm());
 
     return
         gSum
         (
             (mag(Apsi - tmpField) + mag(source - tmpField))(),
-            matrix_.lduMesh_.comm()
+            matrix_.mesh().comm()
         )
       + solverPerformance::small_;
 

src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixUpdateMatrixInterfaces.C

@@ -54,6 +54,8 @@ void Foam::lduMatrix::initMatrixInterfaces
                 (
                     result,
                     add,
+                    mesh().lduAddr(),
+                    interfacei,
                     psiif,
                     coupleCoeffs[interfacei],
                     cmpt,
@@ -81,6 +83,8 @@ void Foam::lduMatrix::initMatrixInterfaces
                 (
                     result,
                     add,
+                    mesh().lduAddr(),
+                    interfacei,
                     psiif,
                     coupleCoeffs[interfacei],
                     cmpt,
@@ -120,6 +124,8 @@ void Foam::lduMatrix::updateMatrixInterfaces
                 (
                     result,
                     add,
+                    mesh().lduAddr(),
+                    interfacei,
                     psiif,
                     coupleCoeffs[interfacei],
                     cmpt,
@@ -149,6 +155,8 @@ void Foam::lduMatrix::updateMatrixInterfaces
                             (
                                 result,
                                 add,
+                                mesh().lduAddr(),
+                                interfacei,
                                 psiif,
                                 coupleCoeffs[interfacei],
                                 cmpt,
@@ -197,6 +205,8 @@ void Foam::lduMatrix::updateMatrixInterfaces
                 (
                     result,
                     add,
+                    mesh().lduAddr(),
+                    interfacei,
                     psiif,
                     coupleCoeffs[interfacei],
                     cmpt,
@@ -210,18 +220,20 @@ void Foam::lduMatrix::updateMatrixInterfaces
         const lduSchedule& patchSchedule = this->patchSchedule();
 
         // Loop over all the "normal" interfaces relating to standard patches
-        forAll(patchSchedule, i)
+        for (const auto& sched : patchSchedule)
         {
-            label interfacei = patchSchedule[i].patch;
+            const label interfacei = sched.patch;
 
             if (interfaces.set(interfacei))
             {
-                if (patchSchedule[i].init)
+                if (sched.init)
                 {
                     interfaces[interfacei].initInterfaceMatrixUpdate
                     (
                         result,
                         add,
+                        mesh().lduAddr(),
+                        interfacei,
                         psiif,
                         coupleCoeffs[interfacei],
                         cmpt,
@@ -234,6 +246,8 @@ void Foam::lduMatrix::updateMatrixInterfaces
                     (
                         result,
                         add,
+                        mesh().lduAddr(),
+                        interfacei,
                         psiif,
                         coupleCoeffs[interfacei],
                         cmpt,
@@ -258,6 +272,8 @@ void Foam::lduMatrix::updateMatrixInterfaces
                 (
                     result,
                     add,
+                    mesh().lduAddr(),
+                    interfacei,
                     psiif,
                     coupleCoeffs[interfacei],
                     cmpt,

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomerateLduAddressing.C

@@ -29,6 +29,7 @@ License
 #include "GAMGAgglomeration.H"
 #include "GAMGInterface.H"
 #include "processorGAMGInterface.H"
+#include "cyclicLduInterface.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -43,7 +44,7 @@ void Foam::GAMGAgglomeration::agglomerateLduAddressing
     const labelUList& upperAddr = fineMeshAddr.upperAddr();
     const labelUList& lowerAddr = fineMeshAddr.lowerAddr();
 
-    label nFineFaces = upperAddr.size();
+    const label nFineFaces = upperAddr.size();
 
     // Get restriction map for current level
     const labelField& restrictMap = restrictAddressing(fineLevelIndex);
@@ -282,15 +283,29 @@ void Foam::GAMGAgglomeration::agglomerateLduAddressing
     const label nReq = Pstream::nRequests();
 
     // Initialise transfer of restrict addressing on the interface
+    // The finest mesh uses patchAddr from the original lduAdressing.
+    // the coarser levels create their own adressing for faceCells
     forAll(fineInterfaces, inti)
     {
         if (fineInterfaces.set(inti))
         {
-            fineInterfaces[inti].initInternalFieldTransfer
-            (
-                Pstream::commsTypes::nonBlocking,
-                restrictMap
-            );
+            if (fineLevelIndex == 0)
+            {
+                fineInterfaces[inti].initInternalFieldTransfer
+                (
+                    Pstream::commsTypes::nonBlocking,
+                    restrictMap,
+                    fineMeshAddr.patchAddr(inti)
+                );
+            }
+            else
+            {
+                fineInterfaces[inti].initInternalFieldTransfer
+                (
+                    Pstream::commsTypes::nonBlocking,
+                    restrictMap
+                );
+            }
         }
     }
 
@@ -320,6 +335,51 @@ void Foam::GAMGAgglomeration::agglomerateLduAddressing
     {
         if (fineInterfaces.set(inti))
         {
+            tmp<labelField> restrictMapInternalField;
+
+            // The finest mesh uses patchAddr from the original lduAdressing.
+            // the coarser levels create thei own adressing for faceCells
+            if (fineLevelIndex == 0)
+            {
+                restrictMapInternalField =
+                    fineInterfaces[inti].interfaceInternalField
+                    (
+                        restrictMap,
+                        fineMeshAddr.patchAddr(inti)
+                    );
+            }
+            else
+            {
+                restrictMapInternalField =
+                    fineInterfaces[inti].interfaceInternalField
+                    (
+                        restrictMap
+                    );
+            }
+
+            tmp<labelField> nbrRestrictMapInternalField =
+                fineInterfaces[inti].internalFieldTransfer
+                (
+                    Pstream::commsTypes::nonBlocking,
+                    restrictMap
+                );
+
+            coarseInterfaces.set
+            (
+                inti,
+                GAMGInterface::New
+                (
+                    inti,
+                    meshLevels_[fineLevelIndex].rawInterfaces(),
+                    fineInterfaces[inti],
+                    restrictMapInternalField(),
+                    nbrRestrictMapInternalField(),
+                    fineLevelIndex,
+                    fineMesh.comm()
+                ).ptr()
+            );
+
+            /* Same as below:
             coarseInterfaces.set
             (
                 inti,
@@ -338,6 +398,7 @@ void Foam::GAMGAgglomeration::agglomerateLduAddressing
                     fineMesh.comm()
                 ).ptr()
             );
+            */
 
             nPatchFaces[inti] = coarseInterfaces[inti].faceCells().size();
             patchFineToCoarse[inti] = refCast<const GAMGInterface>

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C

@@ -295,13 +295,17 @@ const Foam::GAMGAgglomeration& Foam::GAMGAgglomeration::New
     const dictionary& controlDict
 )
 {
-    if
-    (
-        !mesh.thisDb().foundObject<GAMGAgglomeration>
+    const GAMGAgglomeration* agglomPtr =
+        mesh.thisDb().cfindObject<GAMGAgglomeration>
         (
             GAMGAgglomeration::typeName
-        )
-    )
+        );
+
+    if (agglomPtr)
+    {
+        return *agglomPtr;
+    }
+
     {
         const word agglomeratorType
         (
@@ -331,13 +335,6 @@ const Foam::GAMGAgglomeration& Foam::GAMGAgglomeration::New
 
         return store(cstrIter()(mesh, controlDict).ptr());
     }
-    else
-    {
-        return mesh.thisDb().lookupObject<GAMGAgglomeration>
-        (
-            GAMGAgglomeration::typeName
-        );
-    }
 }
 
 
@@ -349,13 +346,17 @@ const Foam::GAMGAgglomeration& Foam::GAMGAgglomeration::New
 {
     const lduMesh& mesh = matrix.mesh();
 
-    if
-    (
-        !mesh.thisDb().foundObject<GAMGAgglomeration>
+    const GAMGAgglomeration* agglomPtr =
+        mesh.thisDb().cfindObject<GAMGAgglomeration>
         (
             GAMGAgglomeration::typeName
-        )
-    )
+        );
+
+    if (agglomPtr)
+    {
+        return *agglomPtr;
+    }
+
     {
         const word agglomeratorType
         (
@@ -385,17 +386,10 @@ const Foam::GAMGAgglomeration& Foam::GAMGAgglomeration::New
             return store(cstrIter()(matrix, controlDict).ptr());
         }
     }
-    else
-    {
-        return mesh.thisDb().lookupObject<GAMGAgglomeration>
-        (
-            GAMGAgglomeration::typeName
-        );
-    }
 }
 
 
-Foam::autoPtr<Foam::GAMGAgglomeration> Foam::GAMGAgglomeration::New
+const Foam::GAMGAgglomeration& Foam::GAMGAgglomeration::New
 (
     const lduMesh& mesh,
     const scalarField& cellVolumes,
@@ -403,40 +397,54 @@ Foam::autoPtr<Foam::GAMGAgglomeration> Foam::GAMGAgglomeration::New
     const dictionary& controlDict
 )
 {
-    const word agglomeratorType
-    (
-        controlDict.getOrDefault<word>("agglomerator", "faceAreaPair")
-    );
 
-    mesh.thisDb().time().libs().open
-    (
-        controlDict,
-        "geometricGAMGAgglomerationLibs",
-        geometryConstructorTablePtr_
-    );
+    const GAMGAgglomeration* agglomPtr =
+        mesh.thisDb().cfindObject<GAMGAgglomeration>
+        (
+            GAMGAgglomeration::typeName
+        );
 
-    auto cstrIter = geometryConstructorTablePtr_->cfind(agglomeratorType);
-
-    if (!cstrIter.found())
+    if (agglomPtr)
     {
-        FatalErrorInFunction
-            << "Unknown GAMGAgglomeration type "
-            << agglomeratorType << ".\n"
-            << "Valid geometric GAMGAgglomeration types :"
-            << geometryConstructorTablePtr_->sortedToc()
-            << exit(FatalError);
+        return *agglomPtr;
     }
 
-    return autoPtr<GAMGAgglomeration>
-    (
-        cstrIter()
+    {
+        const word agglomeratorType
         (
-            mesh,
-            cellVolumes,
-            faceAreas,
-            controlDict
-        )
-    );
+            controlDict.lookupOrDefault<word>("agglomerator", "faceAreaPair")
+        );
+
+        const_cast<Time&>(mesh.thisDb().time()).libs().open
+        (
+            controlDict,
+            "geometricGAMGAgglomerationLibs",
+            geometryConstructorTablePtr_
+        );
+
+        auto cstrIter = geometryConstructorTablePtr_->cfind(agglomeratorType);
+
+        if (!cstrIter.found())
+        {
+            FatalErrorInFunction
+                << "Unknown GAMGAgglomeration type "
+                << agglomeratorType << ".\n"
+                << "Valid geometric GAMGAgglomeration types :"
+                << geometryConstructorTablePtr_->sortedToc()
+                << exit(FatalError);
+        }
+
+        return store
+        (
+            cstrIter()
+            (
+                mesh,
+                cellVolumes,
+                faceAreas,
+                controlDict
+            ).ptr()
+        );
+    }
 }
 
 
@@ -609,8 +617,8 @@ bool Foam::GAMGAgglomeration::checkRestriction
 
         forAll(lower, facei)
         {
-            label own = lower[facei];
-            label nei = upper[facei];
+            const label own = lower[facei];
+            const label nei = upper[facei];
 
             if (restriction[own] == restriction[nei])
             {
@@ -683,7 +691,7 @@ bool Foam::GAMGAgglomeration::checkRestriction
         {
             const label coarseI = restriction[celli];
 
-            label index = coarseToMasters[coarseI].find(master[celli]);
+            const label index = coarseToMasters[coarseI].find(master[celli]);
             newRestrict[celli] = coarseToNewCoarse[coarseI][index];
         }
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.H

@@ -306,7 +306,7 @@ public:
         );
 
         //- Return the selected geometric agglomerator
-        static autoPtr<GAMGAgglomeration> New
+        static const GAMGAgglomeration& New
         (
             const lduMesh& mesh,
             const scalarField& cellVolumes,

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/GAMGInterfaceField/GAMGInterfaceField.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2013 OpenFOAM Foundation
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/cyclicGAMGInterfaceField/cyclicGAMGInterfaceField.C

@@ -107,6 +107,8 @@ void Foam::cyclicGAMGInterfaceField::updateInterfaceMatrix
 (
     solveScalarField& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const solveScalarField& psiInternal,
     const scalarField& coeffs,
     const direction cmpt,
@@ -114,14 +116,20 @@ void Foam::cyclicGAMGInterfaceField::updateInterfaceMatrix
 ) const
 {
     // Get neighbouring field
-    solveScalarField pnf
-    (
-        cyclicInterface_.neighbPatch().interfaceInternalField(psiInternal)
-    );
+
+    const labelList& nbrFaceCells =
+        lduAddr.patchAddr
+        (
+            cyclicInterface_.neighbPatchID()
+        );
+
+    solveScalarField pnf(psiInternal, nbrFaceCells);
 
     transformCoupleField(pnf, cmpt);
 
-    this->addToInternalField(result, !add, coeffs, pnf);
+    const labelList& faceCells = lduAddr.patchAddr(patchId);
+
+    this->addToInternalField(result, !add, faceCells, coeffs, pnf);
 }
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/cyclicGAMGInterfaceField/cyclicGAMGInterfaceField.H

@@ -150,6 +150,8 @@ public:
             (
                 solveScalarField& result,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const solveScalarField& psiInternal,
                 const scalarField& coeffs,
                 const direction cmpt,

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/processorGAMGInterfaceField/processorGAMGInterfaceField.C

@@ -91,6 +91,8 @@ void Foam::processorGAMGInterfaceField::initInterfaceMatrixUpdate
 (
     solveScalarField&,
     const bool,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const solveScalarField& psiInternal,
     const scalarField&,
     const direction,
@@ -142,6 +144,8 @@ void Foam::processorGAMGInterfaceField::updateInterfaceMatrix
 (
     solveScalarField& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const solveScalarField&,
     const scalarField& coeffs,
     const direction cmpt,
@@ -153,6 +157,8 @@ void Foam::processorGAMGInterfaceField::updateInterfaceMatrix
         return;
     }
 
+    const labelUList& faceCells = lduAddr.patchAddr(patchId);
+
     if
     (
         commsType == Pstream::commsTypes::nonBlocking
@@ -178,7 +184,7 @@ void Foam::processorGAMGInterfaceField::updateInterfaceMatrix
         transformCoupleField(scalarReceiveBuf_, cmpt);
 
         // Multiply the field by coefficients and add into the result
-        addToInternalField(result, !add, coeffs, scalarReceiveBuf_);
+        addToInternalField(result, !add, faceCells, coeffs, scalarReceiveBuf_);
     }
     else
     {
@@ -192,7 +198,7 @@ void Foam::processorGAMGInterfaceField::updateInterfaceMatrix
         );
         transformCoupleField(pnf, cmpt);
 
-        addToInternalField(result, !add, coeffs, pnf);
+        addToInternalField(result, !add, faceCells, coeffs, pnf);
     }
 
     const_cast<processorGAMGInterfaceField&>(*this).updatedMatrix() = true;

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/processorGAMGInterfaceField/processorGAMGInterfaceField.H

@@ -142,6 +142,8 @@ public:
             (
                 solveScalarField& result,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const solveScalarField& psiInternal,
                 const scalarField& coeffs,
                 const direction cmpt,
@@ -153,6 +155,8 @@ public:
             (
                 solveScalarField& result,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const solveScalarField& psiInternal,
                 const scalarField& coeffs,
                 const direction cmpt,

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/GAMGInterface/GAMGInterface.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -77,13 +78,23 @@ Foam::tmp<Foam::labelField> Foam::GAMGInterface::interfaceInternalField
 }
 
 
+Foam::tmp<Foam::labelField> Foam::GAMGInterface::interfaceInternalField
+(
+    const labelUList& internalData,
+    const labelUList& faceCells
+) const
+{
+    return interfaceInternalField<label>(internalData, faceCells);
+}
+
+
 Foam::tmp<Foam::scalarField> Foam::GAMGInterface::agglomerateCoeffs
 (
     const scalarField& fineCoeffs
 ) const
 {
-    tmp<scalarField> tcoarseCoeffs(new scalarField(size(), Zero));
-    scalarField& coarseCoeffs = tcoarseCoeffs.ref();
+    auto tcoarseCoeffs = tmp<scalarField>::New(size(), Zero);
+    auto& coarseCoeffs = tcoarseCoeffs.ref();
 
     if (fineCoeffs.size() != faceRestrictAddressing_.size())
     {

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/GAMGInterface/GAMGInterface.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,7 +28,7 @@ Class
     Foam::GAMGInterface
 
 Description
-     Abstract base class for GAMG agglomerated interfaces.
+    Abstract base class for GAMG agglomerated interfaces.
 
 SourceFiles
     GAMGInterface.C
@@ -55,7 +56,6 @@ class GAMGInterface
 :
     public lduInterface
 {
-
 protected:
 
     // Protected data
@@ -73,9 +73,7 @@ protected:
         labelList faceRestrictAddressing_;
 
 
-private:
-
-    // Private Member Functions
+    // Protected Member Functions
 
         //- No copy construct
         GAMGInterface(const GAMGInterface&) = delete;
@@ -245,6 +243,15 @@ public:
                 const UList<Type>& internalData
             ) const;
 
+            //- Return the interface internal field of the given field
+            //- using faceCell mapping
+            template<class Type>
+            tmp<Field<Type>> interfaceInternalField
+            (
+                const UList<Type>& internalData,
+                const labelUList& faceCells
+            ) const;
+
             //- Get the interface internal field of the given field
             template<class Type>
             void interfaceInternalField
@@ -260,6 +267,14 @@ public:
                 const labelUList& internalData
             ) const;
 
+            //- Return the values of the given internal data adjacent to
+            //- the interface as a field using faceCell mapping
+            virtual tmp<labelField> interfaceInternalField
+            (
+                const labelUList& internalData,
+                const labelUList& faceCells
+            ) const;
+
 
         // Agglomeration
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/GAMGInterface/GAMGInterfaceTemplates.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -35,12 +36,30 @@ Foam::tmp<Foam::Field<Type>> Foam::GAMGInterface::interfaceInternalField
     const UList<Type>& iF
 ) const
 {
-    tmp<Field<Type>> tresult(new Field<Type>(size()));
+    auto tresult = tmp<Field<Type>>::New(size());
     interfaceInternalField(iF, tresult.ref());
     return tresult;
 }
 
 
+template<class Type>
+Foam::tmp<Foam::Field<Type>> Foam::GAMGInterface::interfaceInternalField
+(
+    const UList<Type>& iF,
+    const labelUList& faceCells
+) const
+{
+    auto tresult = tmp<Field<Type>>::New(faceCells.size());
+    auto& result = tresult.ref();
+
+    forAll(result, elemi)
+    {
+        result[elemi] = iF[faceCells[elemi]];
+    }
+    return tresult;
+}
+
+
 template<class Type>
 void Foam::GAMGInterface::interfaceInternalField
 (
@@ -48,11 +67,11 @@ void Foam::GAMGInterface::interfaceInternalField
     List<Type>& result
 ) const
 {
-    result.setSize(size());
+    result.resize(size());
 
-    forAll(result, elemI)
+    forAll(result, elemi)
     {
-        result[elemI] = iF[faceCells_[elemI]];
+        result[elemi] = iF[faceCells_[elemi]];
     }
 }
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/cyclicGAMGInterface/cyclicGAMGInterface.C

@@ -166,7 +166,7 @@ Foam::cyclicGAMGInterface::cyclicGAMGInterface
 
 Foam::tmp<Foam::labelField> Foam::cyclicGAMGInterface::internalFieldTransfer
 (
-    const Pstream::commsTypes,
+    const Pstream::commsTypes commsType,
     const labelUList& iF
 ) const
 {

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/cyclicGAMGInterface/cyclicGAMGInterface.H

@@ -162,7 +162,7 @@ public:
         // I/O
 
             //- Write to stream
-            virtual void write(Ostream&) const;
+            virtual void write(Ostream& os) const;
 };
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/processorGAMGInterface/processorGAMGInterface.C

@@ -191,13 +191,24 @@ void Foam::processorGAMGInterface::initInternalFieldTransfer
 }
 
 
+void Foam::processorGAMGInterface::initInternalFieldTransfer
+(
+    const Pstream::commsTypes commsType,
+    const labelUList& iF,
+    const labelUList& faceCells
+) const
+{
+    send(commsType, interfaceInternalField(iF, faceCells)());
+}
+
+
 Foam::tmp<Foam::labelField> Foam::processorGAMGInterface::internalFieldTransfer
 (
     const Pstream::commsTypes commsType,
     const labelUList& iF
 ) const
 {
-    tmp<Field<label>> tfld(receive<label>(commsType, this->size()));
+    tmp<labelField> tfld(receive<label>(commsType, this->size()));
 
     return tfld;
 }

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/processorGAMGInterface/processorGAMGInterface.H

@@ -90,7 +90,7 @@ public:
     // Constructors
 
         //- Construct from fine-level interface,
-        //  local and neighbour restrict addressing
+        //- local and neighbour restrict addressing
         processorGAMGInterface
         (
             const label index,
@@ -140,6 +140,15 @@ public:
                 const labelUList& iF
             ) const;
 
+
+            //- Initialise neighbour field transfer using faceCells mapping
+            virtual void initInternalFieldTransfer
+            (
+                const Pstream::commsTypes commsType,
+                const labelUList& iF,
+                const labelUList& faceCells
+            ) const;
+
             //- Transfer and return internal field adjacent to the interface
             virtual tmp<labelField> internalFieldTransfer
             (
@@ -184,7 +193,7 @@ public:
         // I/O
 
             //- Write to stream
-            virtual void write(Ostream&) const;
+            virtual void write(Ostream& os) const;
 };
 
 

src/OpenFOAM/meshes/lduMesh/lduMesh.H

@@ -89,7 +89,7 @@ public:
             virtual const lduAddressing& lduAddr() const = 0;
 
             //- Return a list of pointers for each patch
-            //  with only those pointing to interfaces being set
+            //- with only those pointing to interfaces being set
             virtual lduInterfacePtrsList interfaces() const = 0;
 
             //- Return communicator used for parallel communication

src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.C

@@ -41,23 +41,24 @@ namespace Foam
     //- Less operator for pairs of \<processor\>\<index\>
     class procLess
     {
-        const labelPairList& lst_;
+        const labelPairList& list_;
 
     public:
 
-        procLess(const labelPairList& lst)
+        procLess(const labelPairList& list)
         :
-            lst_(lst)
+            list_(list)
         {}
 
         bool operator()(const label a, const label b)
         {
-            return lst_[a].first() < lst_[b].first();
+            return list_[a].first() < list_[b].first();
         }
     };
 
 }
 
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::lduPrimitiveMesh::checkUpperTriangular
@@ -85,7 +86,7 @@ void Foam::lduPrimitiveMesh::checkUpperTriangular
         }
     }
 
-    for (label facei=1; facei < l.size(); facei++)
+    for (label facei=1; facei < l.size(); ++facei)
     {
         if (l[facei-1] > l[facei])
         {
@@ -113,6 +114,132 @@ void Foam::lduPrimitiveMesh::checkUpperTriangular
 }
 
 
+Foam::labelListList Foam::lduPrimitiveMesh::globalCellCells
+(
+    const lduMesh& mesh,
+    const globalIndex& globalNumbering
+)
+{
+    const lduAddressing& addr = mesh.lduAddr();
+    lduInterfacePtrsList interfaces = mesh.interfaces();
+
+    const labelList globalIndices
+    (
+        identity
+        (
+            addr.size(),
+            globalNumbering.localStart(UPstream::myProcNo(mesh.comm()))
+        )
+    );
+
+    // Get the interface cells
+    PtrList<labelList> nbrGlobalCells(interfaces.size());
+    {
+        // Initialise transfer of restrict addressing on the interface
+        forAll(interfaces, inti)
+        {
+            if (interfaces.set(inti))
+            {
+                interfaces[inti].initInternalFieldTransfer
+                (
+                    Pstream::commsTypes::nonBlocking,
+                    globalIndices
+                );
+            }
+        }
+
+        if (Pstream::parRun())
+        {
+            Pstream::waitRequests();
+        }
+
+        forAll(interfaces, inti)
+        {
+            if (interfaces.set(inti))
+            {
+                nbrGlobalCells.set
+                (
+                    inti,
+                    new labelList
+                    (
+                        interfaces[inti].internalFieldTransfer
+                        (
+                            Pstream::commsTypes::nonBlocking,
+                            globalIndices
+                        )
+                    )
+                );
+            }
+        }
+    }
+
+
+    // Scan the neighbour list to find out how many times the cell
+    // appears as a neighbour of the face. Done this way to avoid guessing
+    // and resizing list
+    labelList nNbrs(addr.size(), Zero);
+
+    const labelUList& nbr = addr.upperAddr();
+    const labelUList& own = addr.lowerAddr();
+
+    {
+        forAll(nbr, facei)
+        {
+            nNbrs[nbr[facei]]++;
+            nNbrs[own[facei]]++;
+        }
+
+        forAll(interfaces, inti)
+        {
+            if (interfaces.set(inti))
+            {
+                for (const label celli : interfaces[inti].faceCells())
+                {
+                    nNbrs[celli]++;
+                }
+            }
+        }
+    }
+
+
+    // Create cell-cells addressing
+    labelListList cellCells(addr.size());
+
+    forAll(cellCells, celli)
+    {
+        cellCells[celli].setSize(nNbrs[celli], -1);
+    }
+
+    // Reset the list of number of neighbours to zero
+    nNbrs = 0;
+
+    // Scatter the neighbour faces
+    forAll(nbr, facei)
+    {
+        const label c0 = own[facei];
+        const label c1 = nbr[facei];
+
+        cellCells[c0][nNbrs[c0]++] = globalIndices[c1];
+        cellCells[c1][nNbrs[c1]++] = globalIndices[c0];
+    }
+    forAll(interfaces, inti)
+    {
+        if (interfaces.set(inti))
+        {
+            const labelUList& faceCells = interfaces[inti].faceCells();
+
+            forAll(faceCells, facei)
+            {
+                const label c0 = faceCells[facei];
+                cellCells[c0][nNbrs[c0]++] = nbrGlobalCells[inti][facei];
+            }
+        }
+    }
+
+    return cellCells;
+}
+
+
 Foam::label Foam::lduPrimitiveMesh::totalSize
 (
     const PtrList<lduPrimitiveMesh>& meshes
@@ -164,7 +291,7 @@ Foam::labelList Foam::lduPrimitiveMesh::upperTriOrder
     labelList cellToFaces(offsets.last());
     forAll(upper, facei)
     {
-        label celli = lower[facei];
+        const label celli = lower[facei];
         cellToFaces[nNbrs[celli]++] = facei;
     }
 
@@ -172,17 +299,19 @@ Foam::labelList Foam::lduPrimitiveMesh::upperTriOrder
 
     labelList oldToNew(lower.size());
 
-    labelList order;
-    labelList nbr;
+    DynamicList<label> order;
+    DynamicList<label> nbr;
 
     label newFacei = 0;
 
-    for (label celli = 0; celli < nCells; celli++)
+    for (label celli = 0; celli < nCells; ++celli)
     {
-        label startOfCell = offsets[celli];
-        label nNbr = offsets[celli+1] - startOfCell;
+        const label startOfCell = offsets[celli];
+        const label nNbr = offsets[celli+1] - startOfCell;
+
+        nbr.resize(nNbr);
+        order.resize(nNbr);
 
-        nbr.setSize(nNbr);
         forAll(nbr, i)
         {
             nbr[i] = upper[cellToFaces[offsets[celli]+i]];
@@ -216,7 +345,6 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
     comm_(comm)
 {}
 
-
 void Foam::lduPrimitiveMesh::addInterfaces
 (
     lduInterfacePtrsList& interfaces,
@@ -238,6 +366,16 @@ void Foam::lduPrimitiveMesh::addInterfaces
 }
 
 
+Foam::lduPrimitiveMesh::lduPrimitiveMesh
+(
+    const label nCells
+)
+:
+    lduAddressing(nCells),
+    comm_(0)
+{}
+
+
 Foam::lduPrimitiveMesh::lduPrimitiveMesh
 (
     const label nCells,
@@ -251,7 +389,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
     lduAddressing(nCells),
     lowerAddr_(l, true),
     upperAddr_(u, true),
-    primitiveInterfaces_(0),
+    primitiveInterfaces_(),
     patchSchedule_(ps),
     comm_(comm)
 {
@@ -286,10 +424,6 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
 )
 :
     lduAddressing(myMesh.lduAddr().size() + totalSize(otherMeshes)),
-    lowerAddr_(0),
-    upperAddr_(0),
-    interfaces_(0),
-    patchSchedule_(0),
     comm_(comm)
 {
     const label currentComm = myMesh.comm();
@@ -338,7 +472,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
     // Cells get added in order.
     cellOffsets.setSize(nMeshes+1);
     cellOffsets[0] = 0;
-    for (label procMeshI = 0; procMeshI < nMeshes; procMeshI++)
+    for (label procMeshI = 0; procMeshI < nMeshes; ++procMeshI)
     {
         const lduMesh& procMesh = mesh(myMesh, otherMeshes, procMeshI);
 
@@ -351,7 +485,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
     // Faces initially get added in order, sorted later
     labelList internalFaceOffsets(nMeshes+1);
     internalFaceOffsets[0] = 0;
-    for (label procMeshI = 0; procMeshI < nMeshes; procMeshI++)
+    for (label procMeshI = 0; procMeshI < nMeshes; ++procMeshI)
     {
         const lduMesh& procMesh = mesh(myMesh, otherMeshes, procMeshI);
 
@@ -380,7 +514,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
     label nOtherInterfaces = 0;
     labelList nCoupledFaces(nMeshes, Zero);
 
-    for (label procMeshI = 0; procMeshI < nMeshes; procMeshI++)
+    for (label procMeshI = 0; procMeshI < nMeshes; ++procMeshI)
     {
         const lduInterfacePtrsList interfaces =
             mesh(myMesh, otherMeshes, procMeshI).interfaces();
@@ -539,7 +673,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
     faceOffsets.setSize(nMeshes+1);
     faceOffsets[0] = 0;
     faceMap.setSize(nMeshes);
-    for (label procMeshI = 0; procMeshI < nMeshes; procMeshI++)
+    for (label procMeshI = 0; procMeshI < nMeshes; ++procMeshI)
     {
         const lduMesh& procMesh = mesh(myMesh, otherMeshes, procMeshI);
         label nInternal = procMesh.lduAddr().lowerAddr().size();
@@ -566,7 +700,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
     // Old internal faces and resolved coupled interfaces
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    for (label procMeshI = 0; procMeshI < nMeshes; procMeshI++)
+    for (label procMeshI = 0; procMeshI < nMeshes; ++procMeshI)
     {
         const lduMesh& procMesh = mesh(myMesh, otherMeshes, procMeshI);
 
@@ -961,7 +1095,7 @@ void Foam::lduPrimitiveMesh::gather
         otherMeshes.setSize(procIDs.size()-1);
 
         // Slave meshes
-        for (label i = 1; i < procIDs.size(); i++)
+        for (label i = 1; i < procIDs.size(); ++i)
         {
             //Pout<< "on master :"
             //    << " receiving from slave " << procIDs[i] << endl;

src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.H

@@ -28,7 +28,7 @@ Class
     Foam::lduPrimitiveMesh
 
 Description
-    Simplest contrete lduMesh which stores the addressing needed by lduMatrix.
+    Simplest concrete lduMesh that stores the addressing needed by lduMatrix.
 
 SourceFiles
     lduPrimitiveMesh.C
@@ -99,6 +99,9 @@ public:
 
     // Constructors
 
+        //- Construct from number of cells
+        explicit lduPrimitiveMesh(const label nCells);
+
         //- Construct from components but without interfaces. Add interfaces
         //  separately using addInterfaces
         lduPrimitiveMesh
@@ -110,6 +113,7 @@ public:
             bool reuse
         );
 
+
         //- Add interfaces to a mesh constructed without
         void addInterfaces
         (
@@ -186,13 +190,32 @@ public:
                 return *this;
             }
 
+            //- Return non-const ldu addressing
+            virtual lduAddressing& lduAddr()
+            {
+                return *this;
+            }
+
             //- Return a list of pointers for each patch
-            //  with only those pointing to interfaces being set
+            //- with only those pointing to interfaces being set
             virtual lduInterfacePtrsList interfaces() const
             {
                 return interfaces_;
             }
 
+            //- Return a non-const list of pointers for each patch
+            //- with only those pointing to interfaces being set
+            virtual lduInterfacePtrsList& interfaces()
+            {
+                return interfaces_;
+            }
+
+            //- Return a non-const list of primitive interfaces
+            PtrList<const lduInterface>& primitiveInterfaces()
+            {
+                return primitiveInterfaces_;
+            }
+
             //- Return a list of pointers for each patch
             //  with only those pointing to interfaces being set
             //  (reference to cached interfaces)
@@ -219,6 +242,18 @@ public:
                 return upperAddr_;
             }
 
+            //- Return non-const Lower addressing
+            virtual labelList& lowerAddr()
+            {
+                return lowerAddr_;
+            }
+
+            //- Return non-const Upper addressing
+            virtual labelList& upperAddr()
+            {
+                return upperAddr_;
+            }
+
             //- Return patch addressing
             virtual const labelUList& patchAddr(const label i) const
             {
@@ -271,6 +306,13 @@ public:
                 const labelUList& l,
                 const labelUList& u
             );
+
+            //- Calculate global cell-cells
+            static labelListList globalCellCells
+            (
+                const lduMesh& mesh,
+                const globalIndex& globalNumbering
+            );
 };
 
 

src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C

@@ -70,9 +70,9 @@ void Foam::coupledPolyPatch::writeOBJ
     const labelList& pointLabels
 )
 {
-    forAll(pointLabels, i)
+    for (const label pointi : pointLabels)
     {
-        writeOBJ(os, points[pointLabels[i]]);
+        writeOBJ(os, points[pointi]);
     }
 }
 
@@ -108,10 +108,8 @@ void Foam::coupledPolyPatch::writeOBJ
 
     label vertI = 0;
 
-    forAll(faces, i)
+    for (const face& f : faces)
     {
-        const face& f = faces[i];
-
         forAll(f, fp)
         {
             if (foamToObj.insert(f[fp], vertI))
@@ -526,6 +524,18 @@ Foam::coupledPolyPatch::coupledPolyPatch
 {}
 
 
+Foam::coupledPolyPatch::coupledPolyPatch
+(
+    const coupledPolyPatch& pp,
+    const labelList& faceCells
+)
+:
+    polyPatch(pp, faceCells),
+    matchTolerance_(pp.matchTolerance_),
+    transform_(pp.transform_)
+{}
+
+
 Foam::coupledPolyPatch::coupledPolyPatch
 (
     const coupledPolyPatch& pp,

src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -207,11 +208,14 @@ public:
             const word& patchType
         );
 
-        //- Construct as copy, resetting the boundary mesh
-        coupledPolyPatch(const coupledPolyPatch&, const polyBoundaryMesh&);
+        //- Copy construct, resetting the boundary mesh
+        coupledPolyPatch(const coupledPolyPatch&, const polyBoundaryMesh& bm);
+
+        //- Copy construct, resetting the faceCells
+        coupledPolyPatch(const coupledPolyPatch&, const labelList& faceCells);
 
         //- Construct given the original patch and resetting the
-        //  face list and boundary mesh information
+        //- face list and boundary mesh information
         coupledPolyPatch
         (
             const coupledPolyPatch& pp,
@@ -238,6 +242,7 @@ public:
 
     // Member Functions
 
+
         // Access
 
             //- Return true because this patch is coupled
@@ -360,7 +365,7 @@ public:
         );
 
         //- Write the polyPatch data as a dictionary
-        virtual void write(Ostream&) const;
+        virtual void write(Ostream& os) const;
 };
 
 

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C

@@ -734,6 +734,28 @@ Foam::cyclicPolyPatch::cyclicPolyPatch
 }
 
 
+Foam::cyclicPolyPatch::cyclicPolyPatch
+(
+    const cyclicPolyPatch& pp,
+    const label nrbPatchID,
+    const labelList& faceCells
+)
+:
+    coupledPolyPatch(pp, faceCells),
+    neighbPatchName_(pp.neighbPatchName_),
+    coupleGroup_(pp.coupleGroup_),
+    neighbPatchID_(nrbPatchID),
+    rotationAxis_(pp.rotationAxis_),
+    rotationCentre_(pp.rotationCentre_),
+    separationVector_(pp.separationVector_),
+    coupledPointsPtr_(nullptr),
+    coupledEdgesPtr_(nullptr)
+{
+    // Neighbour patch might not be valid yet so no transformation
+    // calculation possible.
+}
+
+
 Foam::cyclicPolyPatch::cyclicPolyPatch
 (
     const cyclicPolyPatch& pp,

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -234,11 +235,19 @@ public:
             const word& patchType
         );
 
-        //- Construct as copy, resetting the boundary mesh
+        //- Copy construct, resetting the boundary mesh
         cyclicPolyPatch(const cyclicPolyPatch&, const polyBoundaryMesh&);
 
+        //- Copy construct, resetting nbrPatchID and faceCells
+        cyclicPolyPatch
+        (
+            const cyclicPolyPatch&,
+            label nbrPatchID,
+            const labelList& faceCells
+        );
+
         //- Construct given the original patch and resetting the
-        //  face list and boundary mesh information
+        //- face list and boundary mesh information
         cyclicPolyPatch
         (
             const cyclicPolyPatch& pp,
@@ -312,6 +321,49 @@ public:
 
     // Member Functions
 
+        // Implicit Functions
+
+             //- Return number of new internal of this polyPatch faces
+            virtual void newInternalProcFaces
+            (
+                label& newFaces,
+                label& newProcFaces
+            ) const
+            {
+                newFaces = this->size();
+                newProcFaces = neighbPatch().size() - this->size();
+            }
+
+            //- Return nbrCells
+            virtual const labelUList& nbrCells() const
+            {
+                return neighbPatch().faceCells();
+            }
+
+            virtual label neighbPolyPatchID() const
+            {
+                return this->neighbPatchID();
+            }
+
+            virtual refPtr<labelListList> mapCollocatedFaces() const
+            {
+                refPtr<labelListList> tMap(new labelListList(this->size()));
+                labelListList& map = tMap.ref();
+                forAll (map, i)
+                {
+                    labelList& subMap = map[i];
+                    subMap.setSize(1);
+                    subMap[0] = i;
+                }
+                return tMap;
+            }
+
+            //- Return implicit master
+            virtual bool masterImplicit() const
+            {
+                return owner();
+            }
+
         //- Neighbour patch name
         const word& neighbPatchName() const;
 

src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C

@@ -252,6 +252,18 @@ Foam::polyPatch::polyPatch(const polyPatch& p)
 {}
 
 
+Foam::polyPatch::polyPatch
+(
+    const polyPatch& p,
+    const labelList& faceCells
+)
+:
+    polyPatch(p)
+{
+    faceCellsPtr_ = new labelList::subList(faceCells, faceCells.size());
+}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::polyPatch::~polyPatch()

src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H

@@ -203,11 +203,11 @@ public:
             const word& patchType
         );
 
-        //- Construct as copy, resetting the boundary mesh
+        //- Copy construct, resetting the boundary mesh
         polyPatch(const polyPatch&, const polyBoundaryMesh&);
 
         //- Construct given the original patch and resetting the
-        //  face list and boundary mesh information
+        //- face list and boundary mesh information
         polyPatch
         (
             const polyPatch& pp,
@@ -227,8 +227,19 @@ public:
             const label newStart
         );
 
-        //- Construct as copy
-        polyPatch(const polyPatch&);
+        //- Copy construct
+        polyPatch(const polyPatch& p);
+
+        //- Copy construct, setting faceCells
+        //  \note The faceCells are held by the caller
+        polyPatch(const polyPatch& p, const labelList& faceCells);
+
+        //- Construct and return a clone, setting faceCells
+        //  \note The faceCells are held by the caller
+        virtual autoPtr<polyPatch> clone(const labelList& faceCells) const
+        {
+            return autoPtr<polyPatch>::New(*this, faceCells);
+        }
 
         //- Construct and return a clone, resetting the boundary mesh
         virtual autoPtr<polyPatch> clone(const polyBoundaryMesh& bm) const
@@ -266,8 +277,7 @@ public:
 
     // Selectors
 
-        //- Return a pointer to a new patch created on freestore from
-        //  components
+        //- Return pointer to a new patch created on freestore from components
         static autoPtr<polyPatch> New
         (
             const word& patchType,
@@ -278,8 +288,7 @@ public:
             const polyBoundaryMesh& bm
         );
 
-        //- Return a pointer to a new patch created on freestore from
-        //  dictionary
+        //- Return pointer to a new patch created on freestore from dictionary
         static autoPtr<polyPatch> New
         (
             const word& name,
@@ -288,8 +297,7 @@ public:
             const polyBoundaryMesh& bm
         );
 
-        //- Return a pointer to a new patch created on freestore from
-        //  dictionary
+        //- Return pointer to a new patch created on freestore from dictionary
         static autoPtr<polyPatch> New
         (
             const word& patchType,
@@ -306,6 +314,51 @@ public:
 
     // Member Functions
 
+
+        // Implicit treatment functions
+
+            //- Return number of new internal of this polyPatch faces
+            virtual void newInternalProcFaces(label&, label&) const
+            {
+                NotImplemented;
+            }
+
+            //- Return nbrCells
+            virtual const labelUList& nbrCells() const
+            {
+                NotImplemented
+                return faceCells();
+            }
+
+            //- Return nbr patchID
+            virtual label neighbPolyPatchID() const
+            {
+                NotImplemented;
+                return -1;
+            }
+
+            //- Return mapped collocated faces
+            virtual refPtr<labelListList> mapCollocatedFaces() const
+            {
+                NotImplemented;
+                refPtr<labelListList> tMap(new labelListList(size()));
+                return tMap;
+            }
+
+            //- Return implicit master
+            virtual bool masterImplicit() const
+            {
+                NotImplemented;
+                return false;
+            }
+
+            //- Return neighbour region name
+            virtual word neighbRegionID() const
+            {
+                return word("none");
+            }
+
+
         //- The offset where this patch starts in the boundary face list
         //  The value is the same as  patch.start() - mesh.nInternalFaces()
         label offset() const;
@@ -364,7 +417,7 @@ public:
 
 
         //- Write the polyPatch data as a dictionary
-        virtual void write(Ostream&) const;
+        virtual void write(Ostream& os) const;
 
 
         // Geometric data; point list required

src/TurbulenceModels/compressible/Make/options

@@ -2,6 +2,7 @@ EXE_INC = \
     -I../turbulenceModels/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/regionCoupled/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \

src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/turbulentTemperatureCoupledBaffleMixed/turbulentTemperatureCoupledBaffleMixedFvPatchScalarField.C

@@ -118,6 +118,23 @@ turbulentTemperatureCoupledBaffleMixedFvPatchScalarField
     kappaLayers_(0),
     contactRes_(0.0)
 {
+    if (!isA<mappedPatchBase>(this->patch().patch()))
+    {
+        FatalErrorInFunction
+            << "' not type '" << mappedPatchBase::typeName << "'"
+            << "\n    for patch " << p.name()
+            << " of field " << internalField().name()
+            << " in file " << internalField().objectPath()
+            << exit(FatalError);
+    }
+
+    WarningInFunction
+        << "This BC has been superseded by "
+        << "compressible::turbulentTemperatureRadCoupledMixed "
+        << "which has more functionalities and it can handle "
+        << "the assemble coupled option for energy. "
+        << endl;
+
     if (dict.readIfPresent("thicknessLayers", thicknessLayers_))
     {
         dict.readEntry("kappaLayers", kappaLayers_);
@@ -349,6 +366,53 @@ void turbulentTemperatureCoupledBaffleMixedFvPatchScalarField::updateCoeffs()
 }
 
 
+void turbulentTemperatureCoupledBaffleMixedFvPatchScalarField::manipulateMatrix
+(
+    fvMatrix<scalar>& m,
+    const label iMatrix,
+    const direction cmpt
+)
+{
+    FatalErrorInFunction
+        << "This BC does not support energy coupling "
+        << "Use compressible::turbulentTemperatureRadCoupledMixed "
+        << "which has more functionalities and it can handle "
+        << "the assemble coupled option for energy. "
+        << abort(FatalError);
+}
+
+
+tmp<Field<scalar>>
+turbulentTemperatureCoupledBaffleMixedFvPatchScalarField::coeffs
+(
+    fvMatrix<scalar>& matrix,
+    const Field<scalar>& coeffs,
+    const label mat
+) const
+{
+    FatalErrorInFunction
+        << "This BC does not support energy coupling "
+        << "Use compressible::turbulentTemperatureRadCoupledMixed "
+        << "which has more functionalities and it can handle "
+        << "the assemble coupled option for energy. "
+        << abort(FatalError);
+    /*
+    const label index(this->patch().index());
+
+    const label nSubFaces(matrix.lduMesh().cellBoundMap()[mat][index].size());
+
+    Field<scalar> mapCoeffs(nSubFaces, Zero);
+
+    label subFaceI = 0;
+    forAll(*this, faceI)
+    {
+
+    }
+    */
+    return tmp<Field<scalar>>(new Field<scalar>());
+}
+
+
 void turbulentTemperatureCoupledBaffleMixedFvPatchScalarField::write
 (
     Ostream& os

src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/turbulentTemperatureCoupledBaffleMixed/turbulentTemperatureCoupledBaffleMixedFvPatchScalarField.H

@@ -117,6 +117,17 @@ class turbulentTemperatureCoupledBaffleMixedFvPatchScalarField
         scalar contactRes_;
 
 
+    // Private member functions
+
+        //- Calculate coefficients for assembly matrix
+        tmp<Field<scalar>> coeffs
+        (
+            fvMatrix<scalar>& matrix,
+            const Field<scalar>&,
+            const label
+        ) const;
+
+
 public:
 
     //- Runtime type information
@@ -193,13 +204,22 @@ public:
         }
 
 
-    // Member functions
+    // Member Functions
 
         //- Update the coefficients associated with the patch field
         virtual void updateCoeffs();
 
+        //- Manipulate matrix
+        virtual void manipulateMatrix
+        (
+            fvMatrix<scalar>& m,
+            const label iMatrix,
+            const direction cmpt
+        );
+
+
         //- Write
-        virtual void write(Ostream&) const;
+        virtual void write(Ostream& os) const;
 };
 
 

src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/turbulentTemperatureRadCoupledMixed/turbulentTemperatureRadCoupledMixedFvPatchScalarField.C

@@ -69,6 +69,7 @@ turbulentTemperatureRadCoupledMixedFvPatchScalarField
     this->refValue() = 0.0;
     this->refGrad() = 0.0;
     this->valueFraction() = 1.0;
+    this->source() = 0.0;
 }
 
 
@@ -152,6 +153,19 @@ turbulentTemperatureRadCoupledMixedFvPatchScalarField
         refGrad() = 0.0;
         valueFraction() = 1.0;
     }
+
+    if (dict.found("useImplicit"))
+    {
+        useImplicit() = dict.get<Switch>("useImplicit");
+    }
+    if (dict.found("source"))
+    {
+        source() = scalarField("source", dict, p.size());
+    }
+    else
+    {
+        source() = 0.0;
+    }
 }
 
 
@@ -320,6 +334,19 @@ void turbulentTemperatureRadCoupledMixedFvPatchScalarField::updateCoeffs()
         refGrad() = (qr + qrNbr)/kappa(Tp);
     }
 
+    source() = Zero;
+
+    // If useImplicit is true we need the source term associated with this BC
+    if (this->useImplicit())
+    {
+        source() =
+            alphaSfDelta()*
+            (
+                valueFraction()*deltaH()
+              + (qr + qrNbr)/beta()
+            );
+    }
+
     mixedFvPatchScalarField::updateCoeffs();
 
     if (debug)
@@ -345,6 +372,138 @@ void turbulentTemperatureRadCoupledMixedFvPatchScalarField::updateCoeffs()
 }
 
 
+void turbulentTemperatureRadCoupledMixedFvPatchScalarField::manipulateMatrix
+(
+    fvMatrix<scalar>& m,
+    const label iMatrix,
+    const direction cmpt
+)
+{
+    FatalErrorInFunction
+        << "This T BC does not support energy coupling "
+        << "It is implemented on he field "
+        << abort(FatalError);
+}
+
+
+tmp<Field<scalar>> turbulentTemperatureRadCoupledMixedFvPatchScalarField::coeffs
+(
+    fvMatrix<scalar>& matrix,
+    const Field<scalar>& coeffs,
+    const label mat
+) const
+{
+    FatalErrorInFunction
+        << "This BC does not support energy coupling "
+        << "Use compressible::turbulentTemperatureRadCoupledMixed "
+        << "which has more functionalities and it can handle "
+        << "the assemble coupled option for energy. "
+        << abort(FatalError);
+
+    return tmp<Field<scalar>>(new Field<scalar>());
+}
+
+
+tmp<scalarField>
+turbulentTemperatureRadCoupledMixedFvPatchScalarField::alphaSfDelta() const
+{
+    return (alpha(*this)*patch().deltaCoeffs()*patch().magSf());
+}
+
+
+tmp<scalarField> turbulentTemperatureRadCoupledMixedFvPatchScalarField::
+beta() const
+{
+    const mappedPatchBase& mpp =
+        refCast<const mappedPatchBase>(patch().patch());
+
+    const label samplePatchi = mpp.samplePolyPatch().index();
+    const polyMesh& nbrMesh = mpp.sampleMesh();
+
+    const fvPatch& nbrPatch =
+        refCast<const fvMesh>(nbrMesh).boundary()[samplePatchi];
+
+    const turbulentTemperatureRadCoupledMixedFvPatchScalarField&
+        nbrField = refCast
+            <const turbulentTemperatureRadCoupledMixedFvPatchScalarField>
+            (
+                nbrPatch.lookupPatchField<volScalarField, scalar>(TnbrName_)
+            );
+
+    // Swap to obtain full local values of neighbour internal field
+    scalarField TcNbr(nbrField.patchInternalField());
+    mpp.distribute(TcNbr);
+
+    scalarField alphaDeltaNbr
+    (
+        nbrField.alpha(TcNbr)*nbrPatch.deltaCoeffs()
+    );
+    mpp.distribute(alphaDeltaNbr);
+
+    scalarField alphaDelta
+    (
+         alpha(*this)*patch().deltaCoeffs()
+    );
+
+    return (alphaDeltaNbr + alphaDelta);
+}
+
+
+tmp<scalarField> turbulentTemperatureRadCoupledMixedFvPatchScalarField::
+deltaH() const
+{
+
+    const mappedPatchBase& mpp =
+        refCast<const mappedPatchBase>(patch().patch());
+
+    const polyMesh& nbrMesh = mpp.sampleMesh();
+
+    const basicThermo* nbrThermo =
+        nbrMesh.cfindObject<basicThermo>(basicThermo::dictName);
+
+    const polyMesh& mesh = patch().boundaryMesh().mesh();
+
+    const basicThermo* localThermo =
+        mesh.cfindObject<basicThermo>(basicThermo::dictName);
+
+
+    if (nbrThermo && localThermo)
+    {
+        const label patchi = patch().index();
+        const scalarField& pp = localThermo->p().boundaryField()[patchi];
+        const scalarField& Tp = *this;
+
+        const mappedPatchBase& mpp =
+            refCast<const mappedPatchBase>(patch().patch());
+
+        const label patchiNrb = mpp.samplePolyPatch().index();
+
+        const scalarField& ppNbr = nbrThermo->p().boundaryField()[patchiNrb];
+        //const scalarField& TpNbr = nbrThermo->T().boundaryField()[patchiNrb];
+
+        // Use this Tp to evaluate he jump as this is updated while doing
+        // updateCoeffs on boundary fields which set T values on patches
+        // then non consistent Tp and Tpnbr could be used from different
+        // updated values (specially when T changes drastically between time
+        // steps/
+        return
+        (
+          -  localThermo->he(pp, Tp, patchi)
+          +  nbrThermo->he(ppNbr, Tp, patchiNrb)
+        );
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Can't find thermos on mapped patch "
+            << " method, but thermo package not available"
+            << exit(FatalError);
+    }
+
+    return tmp<scalarField>::New(patch().size(), Zero);
+}
+
+
 void turbulentTemperatureRadCoupledMixedFvPatchScalarField::write
 (
     Ostream& os

src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/turbulentTemperatureRadCoupledMixed/turbulentTemperatureRadCoupledMixedFvPatchScalarField.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2017 OpenCFD Ltd
+    Copyright (C) 2017-2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -100,12 +100,13 @@ class turbulentTemperatureRadCoupledMixedFvPatchScalarField
     public mixedFvPatchScalarField,
     public temperatureCoupledBase
 {
-    // Private data
+
+    // Private Data
 
         //- Name of field on the neighbour region
         const word TnbrName_;
 
-         //- Name of the radiative heat flux in the neighbour region
+        //- Name of the radiative heat flux in the neighbour region
         const word qrNbrName_;
 
         //- Name of the radiative heat flux in local region
@@ -124,6 +125,27 @@ class turbulentTemperatureRadCoupledMixedFvPatchScalarField
         Switch thermalInertia_;
 
 
+     // Private Functions
+
+            //- Return local alphaSfDelta
+            tmp<scalarField> alphaSfDelta() const;
+
+            //- Return delta enthalpy between regions
+            tmp<scalarField> deltaH() const;
+
+            //- Return the sum of deltaCoeff*alpha from nbr and local
+            tmp<scalarField> beta() const;
+
+
+            //- Calculate coefficients for assembly matrix
+            tmp<Field<scalar>> coeffs
+            (
+                fvMatrix<scalar>& matrix,
+                const Field<scalar>&,
+                const label
+            ) const;
+
+
 public:
 
     //- Runtime type information
@@ -195,13 +217,22 @@ public:
         }
 
 
-    // Member functions
+    // Member Functions
 
         //- Update the coefficients associated with the patch field
         virtual void updateCoeffs();
 
+        //- Manipulate matrix
+        virtual void manipulateMatrix
+        (
+            fvMatrix<scalar>& m,
+            const label iMatrix,
+            const direction cmpt
+        );
+
+
         //- Write
-        virtual void write(Ostream&) const;
+        virtual void write(Ostream& os) const;
 };
 
 

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/porousBafflePressure/porousBafflePressureFvPatchField.C

@@ -130,8 +130,7 @@ void Foam::porousBafflePressureFvPatchField::updateCoeffs()
         return;
     }
 
-    const surfaceScalarField& phi =
-            db().lookupObject<surfaceScalarField>(phiName_);
+    const auto& phi = db().lookupObject<surfaceScalarField>(phiName_);
 
     const fvsPatchField<scalar>& phip =
         patch().patchField<surfaceScalarField, scalar>(phi);

src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -48,6 +49,7 @@ void Foam::cyclicFaPatch::calcTransforms()
 {
     if (size() > 0)
     {
+        // const label sizeby2 = this->size()/2;
         pointField half0Ctrs(size()/2);
         pointField half1Ctrs(size()/2);
         for (label i=0; i<size()/2; ++i)
@@ -153,7 +155,7 @@ void Foam::cyclicFaPatch::makeWeights(scalarField& w) const
     const scalarField& magL = magEdgeLengths();
 
     const scalarField deltas(edgeNormals() & faPatch::delta());
-    label sizeby2 = deltas.size()/2;
+    const label sizeby2 = deltas.size()/2;
 
     scalar maxMatchError = 0;
     label errorEdge = -1;
@@ -206,7 +208,7 @@ void Foam::cyclicFaPatch::makeWeights(scalarField& w) const
 void Foam::cyclicFaPatch::makeDeltaCoeffs(scalarField& dc) const
 {
     const scalarField deltas(edgeNormals() & faPatch::delta());
-    label sizeby2 = deltas.size()/2;
+    const label sizeby2 = deltas.size()/2;
 
     for (label edgei = 0; edgei < sizeby2; ++edgei)
     {
@@ -248,10 +250,10 @@ void Foam::cyclicFaPatch::movePoints(const pointField& p)
 Foam::tmp<Foam::vectorField> Foam::cyclicFaPatch::delta() const
 {
     const vectorField patchD(faPatch::delta());
-    label sizeby2 = patchD.size()/2;
+    const label sizeby2 = patchD.size()/2;
 
-    tmp<vectorField> tpdv(new vectorField(patchD.size()));
-    vectorField& pdv = tpdv.ref();
+    auto tpdv = tmp<vectorField>::New(patchD.size());
+    auto& pdv = tpdv.ref();
 
     // Do the transformation if necessary
     if (parallel())
@@ -290,16 +292,26 @@ Foam::tmp<Foam::labelField> Foam::cyclicFaPatch::interfaceInternalField
 }
 
 
+Foam::tmp<Foam::labelField> Foam::cyclicFaPatch::interfaceInternalField
+(
+    const labelUList& internalData,
+    const labelUList& edgeFaces
+) const
+{
+    return patchInternalField(internalData, edgeFaces);
+}
+
+
 Foam::tmp<Foam::labelField> Foam::cyclicFaPatch::transfer
 (
-    const Pstream::commsTypes,
+    const Pstream::commsTypes commsType,
     const labelUList& interfaceData
 ) const
 {
-    tmp<labelField> tpnf(new labelField(this->size()));
-    labelField& pnf = tpnf.ref();
+    auto tpnf = tmp<labelField>::New(this->size());
+    auto& pnf = tpnf.ref();
 
-    label sizeby2 = this->size()/2;
+    const label sizeby2 = this->size()/2;
 
     for (label edgei=0; edgei<sizeby2; ++edgei)
     {
@@ -317,12 +329,21 @@ Foam::tmp<Foam::labelField> Foam::cyclicFaPatch::internalFieldTransfer
     const labelUList& iF
 ) const
 {
-    const labelUList& edgeCells = this->faceCells();
+    return internalFieldTransfer(commsType, iF, this->faceCells());
+}
 
-    tmp<labelField> tpnf(new labelField(this->size()));
-    labelField& pnf = tpnf.ref();
 
-    label sizeby2 = this->size()/2;
+Foam::tmp<Foam::labelField> Foam::cyclicFaPatch::internalFieldTransfer
+(
+    const Pstream::commsTypes commsType,
+    const labelUList& iF,
+    const labelUList& edgeCells
+) const
+{
+    auto tpnf = tmp<labelField>::New(this->size());
+    auto& pnf = tpnf.ref();
+
+    const label sizeby2 = this->size()/2;
 
     for (label edgei=0; edgei<sizeby2; ++edgei)
     {

src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.H

@@ -167,12 +167,22 @@ public:
         // Interface transfer functions
 
             //- Return the values of the given internal data adjacent to
-            //  the interface as a field
+            //- the interface as a field
             virtual tmp<labelField> interfaceInternalField
             (
                 const labelUList& internalData
             ) const;
 
+
+            //- Return the values of the given internal data adjacent to
+            //- the interface as a field using edgeFace mapping
+            virtual tmp<labelField> interfaceInternalField
+            (
+                const labelUList& internalData,
+                const labelUList& edgeFaces
+            ) const;
+
+
             //- Transfer and return neighbour field
             virtual tmp<labelField> transfer
             (
@@ -186,6 +196,15 @@ public:
                 const Pstream::commsTypes commsType,
                 const labelUList& internalData
             ) const;
+
+
+            //- Return neighbour field using edgeCells mapping
+            virtual tmp<labelField> internalFieldTransfer
+            (
+                const Pstream::commsTypes commsType,
+                const labelUList& internalData,
+                const labelUList& edgeCells
+            ) const;
 };
 
 

src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -77,7 +78,6 @@ void Foam::processorFaPatch::makeNonGlobalPatchPoints() const
     (
         !Pstream::parRun()
      || !boundaryMesh().mesh()().globalData().nGlobalPoints()
-//     || !boundaryMesh().mesh().globalData().nGlobalPoints()
     )
     {
         nonGlobalPatchPointsPtr_ = new labelList(nPoints());
@@ -464,6 +464,16 @@ Foam::tmp<Foam::labelField> Foam::processorFaPatch::interfaceInternalField
 }
 
 
+Foam::tmp<Foam::labelField> Foam::processorFaPatch::interfaceInternalField
+(
+    const labelUList& internalData,
+    const labelUList& edgeFaces
+) const
+{
+    return patchInternalField(internalData, edgeFaces);
+}
+
+
 void Foam::processorFaPatch::initTransfer
 (
     const Pstream::commsTypes commsType,
@@ -504,6 +514,17 @@ Foam::tmp<Foam::labelField> Foam::processorFaPatch::internalFieldTransfer
 }
 
 
+Foam::tmp<Foam::labelField> Foam::processorFaPatch::internalFieldTransfer
+(
+    const Pstream::commsTypes commsType,
+    const labelUList&,
+    const labelUList&
+) const
+{
+    return receive<label>(commsType, this->size());
+}
+
+
 void Foam::processorFaPatch::write(Ostream& os) const
 {
     faPatch::write(os);

src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -264,12 +265,20 @@ public:
         // Interface transfer functions
 
             //- Return the values of the given internal data adjacent to
-            //  the interface as a field
+            //- the interface as a field
             virtual tmp<labelField> interfaceInternalField
             (
                 const labelUList& internalData
             ) const;
 
+            //- Return the values of the given internal data adjacent to
+            //- the interface as a field using edgeFaces
+            virtual tmp<labelField> interfaceInternalField
+            (
+                const labelUList& internalData,
+                const labelUList& edgeFaces
+            ) const;
+
             //- Initialise interface data transfer
             virtual void initTransfer
             (
@@ -298,8 +307,16 @@ public:
                 const labelUList& internalData
             ) const;
 
+            //- Return neighbour field using mapping
+            virtual tmp<labelField> internalFieldTransfer
+            (
+                const Pstream::commsTypes commsType,
+                const labelUList& internalData,
+                const labelUList& edgeCells
+            ) const;
+
             //- Write the patch data as a dictionary
-            virtual void write(Ostream&) const;
+            virtual void write(Ostream& os) const;
 };
 
 

src/finiteArea/faMesh/faPatches/faPatch/faPatch.H

@@ -357,6 +357,15 @@ public:
             template<class Type>
             tmp<Field<Type>> patchInternalField(const UList<Type>&) const;
 
+            //- Return given internal field next to patch as patch field
+            //  providing addressing
+            template<class Type>
+            tmp<Foam::Field<Type>> patchInternalField
+            (
+                const UList<Type>& f,
+                const labelUList& edgeFaces
+            ) const;
+
             //- Return the corresponding patchField of the named field
             template<class GeometricField, class Type>
             const typename GeometricField::Patch& patchField

src/finiteArea/faMesh/faPatches/faPatch/faPatchTemplates.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -35,10 +36,19 @@ Foam::tmp<Foam::Field<Type>> Foam::faPatch::patchInternalField
     const UList<Type>& f
 ) const
 {
-    tmp<Field<Type>> tpif (new Field<Type>(size()));
-    Field<Type>& pif = tpif.ref();
+    return patchInternalField(f, this->edgeFaces());
+}
 
-    const labelUList& edgeFaces = this->edgeFaces();
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>> Foam::faPatch::patchInternalField
+(
+    const UList<Type>& f,
+    const labelUList& edgeFaces
+) const
+{
+    auto tpif = tmp<Field<Type>>::New(size());
+    auto& pif = tpif.ref();
 
     forAll(pif, facei)
     {

src/finiteArea/fields/faPatchFields/basic/coupled/coupledFaPatchField.H

@@ -193,6 +193,8 @@ public:
             (
                 solveScalarField& result,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const solveScalarField& psiInternal,
                 const scalarField& coeffs,
                 const direction,
@@ -204,6 +206,8 @@ public:
             (
                 Field<Type>& result,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const Field<Type>&,
                 const scalarField& coeffs,
                 const Pstream::commsTypes commsType

src/finiteArea/fields/faPatchFields/constraint/cyclic/cyclicFaPatchField.C

@@ -129,7 +129,7 @@ Foam::cyclicFaPatchField<Type>::patchNeighbourField() const
     tmp<Field<Type>> tpnf(new Field<Type>(this->size()));
     Field<Type>& pnf = tpnf.ref();
 
-    label sizeby2 = this->size()/2;
+    const label sizeby2 = this->size()/2;
 
     if (doTransform())
     {
@@ -164,6 +164,8 @@ void Foam::cyclicFaPatchField<Type>::updateInterfaceMatrix
 (
     solveScalarField& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const solveScalarField& psiInternal,
     const scalarField& coeffs,
     const direction cmpt,
@@ -172,7 +174,7 @@ void Foam::cyclicFaPatchField<Type>::updateInterfaceMatrix
 {
     solveScalarField pnf(this->size());
 
-    label sizeby2 = this->size()/2;
+    const label sizeby2 = this->size()/2;
     const labelUList& faceCells = cyclicPatch_.faceCells();
 
     for (label facei = 0; facei < sizeby2; ++facei)
@@ -207,6 +209,8 @@ void Foam::cyclicFaPatchField<Type>::updateInterfaceMatrix
 (
     Field<Type>& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const Field<Type>& psiInternal,
     const scalarField& coeffs,
     const Pstream::commsTypes commsType
@@ -214,7 +218,7 @@ void Foam::cyclicFaPatchField<Type>::updateInterfaceMatrix
 {
     Field<Type> pnf(this->size());
 
-    label sizeby2 = this->size()/2;
+    const label sizeby2 = this->size()/2;
     const labelUList& faceCells = cyclicPatch_.faceCells();
 
     for (label facei = 0; facei < sizeby2; ++facei)

src/finiteArea/fields/faPatchFields/constraint/cyclic/cyclicFaPatchField.H

@@ -179,6 +179,8 @@ public:
             (
                 solveScalarField& result,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const solveScalarField& psiInternal,
                 const scalarField& coeffs,
                 const direction cmpt,
@@ -190,6 +192,8 @@ public:
             (
                 Field<Type>& result,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const Field<Type>&,
                 const scalarField& coeffs,
                 const Pstream::commsTypes commsType

src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchField.C

@@ -192,6 +192,8 @@ void Foam::processorFaPatchField<Type>::initInterfaceMatrixUpdate
 (
     solveScalarField& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const solveScalarField& psiInternal,
     const scalarField& coeffs,
     const direction,
@@ -211,6 +213,8 @@ void Foam::processorFaPatchField<Type>::updateInterfaceMatrix
 (
     solveScalarField& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const solveScalarField&,
     const scalarField& coeffs,
     const direction cmpt,
@@ -251,6 +255,8 @@ void Foam::processorFaPatchField<Type>::initInterfaceMatrixUpdate
 (
     Field<Type>& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const Field<Type>& psiInternal,
     const scalarField& coeffs,
     const Pstream::commsTypes commsType
@@ -269,6 +275,8 @@ void Foam::processorFaPatchField<Type>::updateInterfaceMatrix
 (
     Field<Type>& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const Field<Type>&,
     const scalarField& coeffs,
     const Pstream::commsTypes commsType

src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchField.H

@@ -193,6 +193,8 @@ public:
             (
                 solveScalarField& result,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const solveScalarField& psiInternal,
                 const scalarField& coeffs,
                 const direction cmpt,
@@ -204,6 +206,8 @@ public:
             (
                 solveScalarField& result,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const solveScalarField& psiInternal,
                 const scalarField& coeffs,
                 const direction cmpt,
@@ -215,6 +219,8 @@ public:
             (
                 Field<Type>& result,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const Field<Type>&,
                 const scalarField& coeffs,
                 const Pstream::commsTypes commsType
@@ -225,6 +231,8 @@ public:
             (
                 Field<Type>& result,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const Field<Type>&,
                 const scalarField& coeffs,
                 const Pstream::commsTypes commsType

src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchScalarField.C

@@ -44,6 +44,8 @@ void Foam::processorFaPatchField<Foam::scalar>::initInterfaceMatrixUpdate
 (
     solveScalarField& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const solveScalarField& psiInternal,
     const scalarField& coeffs,
     const direction,
@@ -63,6 +65,8 @@ void Foam::processorFaPatchField<Foam::scalar>::updateInterfaceMatrix
 (
     solveScalarField& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const solveScalarField&,
     const scalarField& coeffs,
     const direction,

src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchScalarField.H

@@ -51,6 +51,8 @@ void processorFaPatchField<scalar>::initInterfaceMatrixUpdate
 (
     solveScalarField& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const solveScalarField&,
     const scalarField& coeffs,
     const direction,
@@ -63,6 +65,8 @@ void processorFaPatchField<scalar>::updateInterfaceMatrix
 (
     solveScalarField& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const solveScalarField&,
     const scalarField& coeffs,
     const direction,

src/finiteVolume/Make/files

@@ -142,6 +142,7 @@ $(constraintFvPatchFields)/jumpCyclicAMI/jumpCyclicAMIFvPatchFields.C
 $(constraintFvPatchFields)/nonuniformTransformCyclic/nonuniformTransformCyclicFvPatchFields.C
 $(constraintFvPatchFields)/processor/processorFvPatchFields.C
 $(constraintFvPatchFields)/processorCyclic/processorCyclicFvPatchFields.C
+$(constraintFvPatchFields)/calculatedProcessor/calculatedProcessorFvPatchFields.C
 $(constraintFvPatchFields)/symmetryPlane/symmetryPlaneFvPatchFields.C
 $(constraintFvPatchFields)/symmetryPlane/symmetryPlaneFvPatchScalarField.C
 $(constraintFvPatchFields)/symmetry/symmetryFvPatchFields.C
@@ -516,4 +517,9 @@ $(fvOptions)/fvOptionIO.C
 $(fvOptions)/fvOptionList.C
 $(fvOptions)/fvOptions.C
 
+lduPrimitiveMeshAssembly/AssemblyFvPatches.C
+lduPrimitiveMeshAssembly/lduPrimitiveMeshAssembly.C
+lduPrimitiveMeshAssembly/assemblyFaceAreaPairGAMGAgglomeration/assemblyFaceAreaPairGAMGAgglomeration.C
+
+
 LIB = $(FOAM_LIBBIN)/libfiniteVolume

src/finiteVolume/Make/options

@@ -9,3 +9,4 @@ LIB_LIBS = \
     -lfileFormats \
     -lsurfMesh \
     -lmeshTools
+

src/finiteVolume/fields/fvPatchFields/basic/coupled/coupledFvPatchField.H

@@ -211,6 +211,8 @@ public:
             (
                 solveScalarField& result,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const solveScalarField& psiInternal,
                 const scalarField& coeffs,
                 const direction,
@@ -222,6 +224,8 @@ public:
             (
                 Field<Type>&,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const Field<Type>&,
                 const scalarField&,
                 const Pstream::commsTypes commsType

src/finiteVolume/fields/fvPatchFields/basic/mixed/mixedFvPatchField.C

@@ -39,7 +39,8 @@ Foam::mixedFvPatchField<Type>::mixedFvPatchField
     fvPatchField<Type>(p, iF),
     refValue_(p.size()),
     refGrad_(p.size()),
-    valueFraction_(p.size())
+    valueFraction_(p.size()),
+    source_(p.size(), Zero)
 {}
 
 
@@ -54,7 +55,8 @@ Foam::mixedFvPatchField<Type>::mixedFvPatchField
     fvPatchField<Type>(p, iF, dict, false),
     refValue_("refValue", dict, p.size()),
     refGrad_("refGradient", dict, p.size()),
-    valueFraction_("valueFraction", dict, p.size())
+    valueFraction_("valueFraction", dict, p.size()),
+    source_(p.size(), Zero)
 {
     // Could also check/clip fraction to 0-1 range
     evaluate();
@@ -73,7 +75,8 @@ Foam::mixedFvPatchField<Type>::mixedFvPatchField
     fvPatchField<Type>(ptf, p, iF, mapper),
     refValue_(ptf.refValue_, mapper),
     refGrad_(ptf.refGrad_, mapper),
-    valueFraction_(ptf.valueFraction_, mapper)
+    valueFraction_(ptf.valueFraction_, mapper),
+    source_(ptf.source_, mapper)
 {
     if (notNull(iF) && mapper.hasUnmapped())
     {
@@ -96,7 +99,8 @@ Foam::mixedFvPatchField<Type>::mixedFvPatchField
     fvPatchField<Type>(ptf),
     refValue_(ptf.refValue_),
     refGrad_(ptf.refGrad_),
-    valueFraction_(ptf.valueFraction_)
+    valueFraction_(ptf.valueFraction_),
+    source_(ptf.source_)
 {}
 
 
@@ -110,7 +114,8 @@ Foam::mixedFvPatchField<Type>::mixedFvPatchField
     fvPatchField<Type>(ptf, iF),
     refValue_(ptf.refValue_),
     refGrad_(ptf.refGrad_),
-    valueFraction_(ptf.valueFraction_)
+    valueFraction_(ptf.valueFraction_),
+    source_(ptf.source_)
 {}
 
 
@@ -126,6 +131,7 @@ void Foam::mixedFvPatchField<Type>::autoMap
     refValue_.autoMap(m);
     refGrad_.autoMap(m);
     valueFraction_.autoMap(m);
+    source_.autoMap(m);
 }
 
 
@@ -144,12 +150,14 @@ void Foam::mixedFvPatchField<Type>::rmap
     refValue_.rmap(mptf.refValue_, addr);
     refGrad_.rmap(mptf.refGrad_, addr);
     valueFraction_.rmap(mptf.valueFraction_, addr);
+    source_.rmap(mptf.source_, addr);
 }
 
 
 template<class Type>
 void Foam::mixedFvPatchField<Type>::evaluate(const Pstream::commsTypes)
 {
+
     if (!this->updated())
     {
         this->updateCoeffs();
@@ -158,9 +166,8 @@ void Foam::mixedFvPatchField<Type>::evaluate(const Pstream::commsTypes)
     Field<Type>::operator=
     (
         valueFraction_*refValue_
-      +
-        (1.0 - valueFraction_)*
-        (
+      + (1.0 - valueFraction_)
+       *(
             this->patchInternalField()
           + refGrad_/this->patch().deltaCoeffs()
         )
@@ -190,7 +197,6 @@ Foam::mixedFvPatchField<Type>::valueInternalCoeffs
 ) const
 {
     return Type(pTraits<Type>::one)*(1.0 - valueFraction_);
-
 }
 
 
@@ -202,8 +208,8 @@ Foam::mixedFvPatchField<Type>::valueBoundaryCoeffs
 ) const
 {
     return
-         valueFraction_*refValue_
-       + (1.0 - valueFraction_)*refGrad_/this->patch().deltaCoeffs();
+        valueFraction_*refValue_
+      + (1.0 - valueFraction_)*refGrad_/this->patch().deltaCoeffs();
 }
 
 
@@ -232,6 +238,7 @@ void Foam::mixedFvPatchField<Type>::write(Ostream& os) const
     refValue_.writeEntry("refValue", os);
     refGrad_.writeEntry("refGradient", os);
     valueFraction_.writeEntry("valueFraction", os);
+    source_.writeEntry("source", os);
     this->writeEntry("value", os);
 }
 

src/finiteVolume/fields/fvPatchFields/basic/mixed/mixedFvPatchField.H

@@ -100,6 +100,9 @@ class mixedFvPatchField
         //- Fraction (0-1) of value used for boundary condition
         scalarField valueFraction_;
 
+        //- Source field
+        Field<Type> source_;
+
 
 public:
 
@@ -219,6 +222,16 @@ public:
                 return valueFraction_;
             }
 
+            virtual Field<Type>& source()
+            {
+                return source_;
+            }
+
+            virtual const Field<Type>& source() const
+            {
+                return source_;
+            }
+
 
         // Mapping functions
 

src/finiteVolume/fields/fvPatchFields/constraint/calculatedProcessor/calculatedProcessorFvPatchField.C

@@ -37,7 +37,6 @@ Foam::calculatedProcessorFvPatchField<Type>::calculatedProcessorFvPatchField
     const DimensionedField<Type, volMesh>& iF
 )
 :
-    //lduInterfaceField(interface),
     coupledFvPatchField<Type>(p, iF),
     procInterface_(refCast<const lduPrimitiveProcessorInterface>(interface)),
     sendBuf_(interface.faceCells().size()),
@@ -55,7 +54,6 @@ Foam::calculatedProcessorFvPatchField<Type>::calculatedProcessorFvPatchField
     const calculatedProcessorFvPatchField<Type>& ptf
 )
 :
-    //lduInterfaceField(ptf.procInterface_),
     coupledFvPatchField<Type>(ptf),
     procInterface_(ptf.procInterface_),
     sendBuf_(procInterface_.faceCells().size()),
@@ -75,7 +73,6 @@ Foam::calculatedProcessorFvPatchField<Type>::calculatedProcessorFvPatchField
     const DimensionedField<Type, volMesh>& iF
 )
 :
-    //lduInterfaceField(ptf.procInterface_),
     coupledFvPatchField<Type>(ptf, iF),
     procInterface_(ptf.procInterface_),
     sendBuf_(procInterface_.faceCells().size()),
@@ -98,9 +95,7 @@ bool Foam::calculatedProcessorFvPatchField<Type>::ready() const
      && this->outstandingSendRequest_ < Pstream::nRequests()
     )
     {
-        bool finished =
-            UPstream::finishedRequest(this->outstandingSendRequest_);
-        if (!finished)
+        if (!UPstream::finishedRequest(this->outstandingSendRequest_))
         {
             return false;
         }
@@ -113,9 +108,7 @@ bool Foam::calculatedProcessorFvPatchField<Type>::ready() const
      && this->outstandingRecvRequest_ < Pstream::nRequests()
     )
     {
-        bool finished =
-            UPstream::finishedRequest(this->outstandingRecvRequest_);
-        if (!finished)
+        if (!UPstream::finishedRequest(this->outstandingRecvRequest_))
         {
             return false;
         }
@@ -130,7 +123,7 @@ template<class Type>
 Foam::tmp<Foam::Field<Type>>
 Foam::calculatedProcessorFvPatchField<Type>::patchNeighbourField() const
 {
-    if (debug && !this->ready())
+    if (!this->ready())
     {
         FatalErrorInFunction
             << "On patch of size " << procInterface_.faceCells().size()
@@ -217,6 +210,8 @@ void Foam::calculatedProcessorFvPatchField<Type>::initInterfaceMatrixUpdate
 (
     solveScalarField& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const solveScalarField& psiInternal,
     const scalarField& coeffs,
     const direction cmpt,
@@ -224,14 +219,16 @@ void Foam::calculatedProcessorFvPatchField<Type>::initInterfaceMatrixUpdate
 ) const
 {
     // Bypass patchInternalField since uses fvPatch addressing
-    const labelList& fc = procInterface_.faceCells();
+
+    const labelList& fc = lduAddr.patchAddr(patchId);
+
     scalarSendBuf_.setSize(fc.size());
     forAll(fc, i)
     {
         scalarSendBuf_[i] = psiInternal[fc[i]];
     }
 
-    if (debug && !this->ready())
+    if (!this->ready())
     {
         FatalErrorInFunction
             << "On patch " //<< interface_.name()
@@ -240,8 +237,10 @@ void Foam::calculatedProcessorFvPatchField<Type>::initInterfaceMatrixUpdate
     }
 
 
+
     scalarReceiveBuf_.setSize(scalarSendBuf_.size());
     outstandingRecvRequest_ = UPstream::nRequests();
+
     UIPstream::read
     (
         Pstream::commsTypes::nonBlocking,
@@ -253,6 +252,7 @@ void Foam::calculatedProcessorFvPatchField<Type>::initInterfaceMatrixUpdate
     );
 
     outstandingSendRequest_ = UPstream::nRequests();
+
     UOPstream::write
     (
         Pstream::commsTypes::nonBlocking,
@@ -303,6 +303,8 @@ void Foam::calculatedProcessorFvPatchField<Type>::updateInterfaceMatrix
 (
     solveScalarField& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const solveScalarField& psiInternal,
     const scalarField& coeffs,
     const direction cmpt,
@@ -314,6 +316,7 @@ void Foam::calculatedProcessorFvPatchField<Type>::updateInterfaceMatrix
         return;
     }
 
+
     if
     (
         outstandingRecvRequest_ >= 0
@@ -337,13 +340,4 @@ void Foam::calculatedProcessorFvPatchField<Type>::updateInterfaceMatrix
 }
 
 
-//template<class Type>
-//void Foam::calculatedProcessorFvPatchField<Type>::write(Ostream& os) const
-//{
-//    //zeroGradientFvPatchField<Type>::write(os);
-//    fvPatchField<Type>::write(os);
-//    this->writeEntry("value", os);
-//}
-
-
 // ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/constraint/calculatedProcessor/calculatedProcessorFvPatchField.H

@@ -30,7 +30,7 @@ Group
     grpGenericBoundaryConditions
 
 Description
-    processorFvPatchField type bypassing fvPatch
+    A processorFvPatchField type bypassing fvPatch
 
     Used to temporarily add updateInterfaceMatrix capabilities to a matrix
     during overset solving. Supplies:
@@ -68,7 +68,7 @@ class calculatedProcessorFvPatchField
 {
 protected:
 
-    // Private data
+    // Protected data
 
          //- Local reference cast into the interface
          const lduPrimitiveProcessorInterface& procInterface_;
@@ -94,8 +94,7 @@ protected:
             mutable label outstandingRecvRequest_;
 
 
-
-    // Private Member Functions
+    // Protected Member Functions
 
         void addToInternalField
         (
@@ -234,6 +233,8 @@ public:
             (
                 solveScalarField& result,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const solveScalarField& psiInternal,
                 const scalarField& coeffs,
                 const direction cmpt,
@@ -245,6 +246,8 @@ public:
             (
                 solveScalarField& result,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const solveScalarField& psiInternal,
                 const scalarField& coeffs,
                 const direction cmpt,
@@ -256,6 +259,8 @@ public:
             (
                 Field<Type>& result,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const Field<Type>& psiInternal,
                 const scalarField& coeffs,
                 const Pstream::commsTypes commsType
@@ -269,6 +274,8 @@ public:
             (
                 Field<Type>& result,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const Field<Type>& psiInternal,
                 const scalarField& coeffs,
                 const Pstream::commsTypes commsType
@@ -278,7 +285,7 @@ public:
             }
 
         //- Write
-        //virtual void write(Ostream&) const;
+        // virtual void write(Ostream& os) const;
 };
 
 

src/finiteVolume/fields/fvPatchFields/constraint/calculatedProcessor/calculatedProcessorFvPatchFields.C

@@ -37,15 +37,31 @@ namespace Foam
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-defineNamedTemplateTypeNameAndDebug(calculatedProcessorFvPatchScalarField, 0);
-defineNamedTemplateTypeNameAndDebug(calculatedProcessorFvPatchVectorField, 0);
+defineNamedTemplateTypeNameAndDebug
+(
+    calculatedProcessorFvPatchScalarField,
+    0
+);
+defineNamedTemplateTypeNameAndDebug
+(
+    calculatedProcessorFvPatchVectorField,
+    0
+);
 defineNamedTemplateTypeNameAndDebug
 (
     calculatedProcessorFvPatchSphericalTensorField,
     0
 );
-defineNamedTemplateTypeNameAndDebug(calculatedProcessorFvPatchSymmTensorField, 0);
-defineNamedTemplateTypeNameAndDebug(calculatedProcessorFvPatchTensorField, 0);
+defineNamedTemplateTypeNameAndDebug
+(
+    calculatedProcessorFvPatchSymmTensorField,
+    0
+);
+defineNamedTemplateTypeNameAndDebug
+(
+    calculatedProcessorFvPatchTensorField,
+    0
+);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.C

@@ -26,8 +26,10 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+#include "fvMatrix.H"
 #include "cyclicFvPatchField.H"
 #include "transformField.H"
+#include "volFields.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -174,27 +176,36 @@ Foam::cyclicFvPatchField<Type>::neighbourPatchField() const
 }
 
 
+
 template<class Type>
 void Foam::cyclicFvPatchField<Type>::updateInterfaceMatrix
 (
     solveScalarField& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const solveScalarField& psiInternal,
     const scalarField& coeffs,
     const direction cmpt,
-    const Pstream::commsTypes
+    const Pstream::commsTypes commsType
 ) const
 {
     const labelUList& nbrFaceCells =
-        cyclicPatch().cyclicPatch().neighbPatch().faceCells();
+        lduAddr.patchAddr
+        (
+            this->cyclicPatch().neighbPatchID()
+        );
 
     solveScalarField pnf(psiInternal, nbrFaceCells);
 
     // Transform according to the transformation tensors
     transformCoupleField(pnf, cmpt);
 
+
+    const labelUList& faceCells = lduAddr.patchAddr(patchId);
+
     // Multiply the field by coefficients and add into the result
-    this->addToInternalField(result, !add, coeffs, pnf);
+    this->addToInternalField(result, !add, faceCells, coeffs, pnf);
 }
 
 
@@ -203,21 +214,28 @@ void Foam::cyclicFvPatchField<Type>::updateInterfaceMatrix
 (
     Field<Type>& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const Field<Type>& psiInternal,
     const scalarField& coeffs,
     const Pstream::commsTypes
 ) const
 {
-    const labelUList& nbrFaceCells =
-        cyclicPatch().cyclicPatch().neighbPatch().faceCells();
+    const labelList& nbrFaceCells =
+        lduAddr.patchAddr
+        (
+            this->cyclicPatch().neighbPatchID()
+        );
 
     Field<Type> pnf(psiInternal, nbrFaceCells);
 
     // Transform according to the transformation tensors
     transformCoupleField(pnf);
 
+    const labelUList& faceCells = lduAddr.patchAddr(patchId);
+
     // Multiply the field by coefficients and add into the result
-    this->addToInternalField(result, !add, coeffs, pnf);
+    this->addToInternalField(result, !add, faceCells, coeffs, pnf);
 }
 
 
@@ -228,4 +246,81 @@ void Foam::cyclicFvPatchField<Type>::write(Ostream& os) const
 }
 
 
+template<class Type>
+void Foam::cyclicFvPatchField<Type>::manipulateMatrix
+(
+    fvMatrix<Type>& matrix,
+    const label mat,
+    const direction cmpt
+)
+{
+    if (this->cyclicPatch().owner())
+    {
+        label index = this->patch().index();
+
+        const label globalPatchID =
+            matrix.lduMeshAssembly().patchLocalToGlobalMap()[mat][index];
+
+        const Field<scalar> intCoeffsCmpt
+        (
+            matrix.internalCoeffs()[globalPatchID].component(cmpt)
+        );
+
+        const Field<scalar> boundCoeffsCmpt
+        (
+            matrix.boundaryCoeffs()[globalPatchID].component(cmpt)
+        );
+
+        const labelUList& u = matrix.lduAddr().upperAddr();
+        const labelUList& l = matrix.lduAddr().lowerAddr();
+
+        const labelList& faceMap =
+            matrix.lduMeshAssembly().faceBoundMap()[mat][index];
+
+        forAll (faceMap, faceI)
+        {
+            label globalFaceI = faceMap[faceI];
+
+            const scalar boundCorr = -boundCoeffsCmpt[faceI];
+            const scalar intCorr = -intCoeffsCmpt[faceI];
+
+            matrix.upper()[globalFaceI] += boundCorr;
+            matrix.diag()[u[globalFaceI]] -= boundCorr;
+            matrix.diag()[l[globalFaceI]] -= intCorr;
+
+            if (matrix.asymmetric())
+            {
+                matrix.lower()[globalFaceI] += intCorr;
+            }
+        }
+
+        if (matrix.psi(mat).mesh().fluxRequired(this->internalField().name()))
+        {
+            matrix.internalCoeffs().set
+            (
+                globalPatchID, intCoeffsCmpt*pTraits<Type>::one
+            );
+            matrix.boundaryCoeffs().set
+            (
+                globalPatchID, boundCoeffsCmpt*pTraits<Type>::one
+            );
+
+            const label nbrPathID = this->cyclicPatch().neighbPatchID();
+
+            const label nbrGlobalPatchID =
+                matrix.lduMeshAssembly().patchLocalToGlobalMap()[mat][nbrPathID];
+
+            matrix.internalCoeffs().set
+            (
+                nbrGlobalPatchID, intCoeffsCmpt*pTraits<Type>::one
+
+            );
+            matrix.boundaryCoeffs().set
+            (
+                nbrGlobalPatchID, boundCoeffsCmpt*pTraits<Type>::one
+            );
+        }
+    }
+}
+
 // ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.H

@@ -93,6 +93,7 @@ class cyclicFvPatchField
 
 public:
 
+
     //- Runtime type information
     TypeName(cyclicFvPatch::typeName_());
 
@@ -183,6 +184,8 @@ public:
             (
                 solveScalarField& result,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const solveScalarField& psiInternal,
                 const scalarField& coeffs,
                 const direction cmpt,
@@ -194,12 +197,24 @@ public:
             (
                 Field<Type>& result,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const Field<Type>& psiInternal,
                 const scalarField& coeffs,
                 const Pstream::commsTypes commsType
             ) const;
 
 
+            //- Manipulate matrix
+            virtual void manipulateMatrix
+            (
+                fvMatrix<Type>& m,
+                const label iMatrix,
+                const direction cmp
+            );
+
+
+
         // Cyclic coupled interface functions
 
             //- Does the patch field perform the transformation

src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.C

@@ -26,6 +26,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+#include "fvMatrix.H"
 #include "cyclicACMIFvPatchField.H"
 #include "transformField.H"
 
@@ -160,7 +161,13 @@ Foam::tmp<Foam::Field<Type>>
 Foam::cyclicACMIFvPatchField<Type>::patchNeighbourField() const
 {
     const Field<Type>& iField = this->primitiveField();
-    const cyclicACMIPolyPatch& cpp = cyclicACMIPatch_.cyclicACMIPatch();
+    //const cyclicACMIPolyPatch& cpp = cyclicACMIPatch_.cyclicACMIPatch();
+
+    // By pass polyPatch to get nbrId. Instead use cyclicAMIFvPatch virtual
+    // neighbPatch()
+    const cyclicACMIFvPatch& neighbPatch = cyclicACMIPatch_.neighbPatch();
+    const labelUList& nbrFaceCells = neighbPatch.faceCells();
+
     tmp<Field<Type>> tpnf
     (
         cyclicACMIPatch_.interpolate
@@ -168,7 +175,8 @@ Foam::cyclicACMIFvPatchField<Type>::patchNeighbourField() const
             Field<Type>
             (
                 iField,
-                cpp.neighbPatch().faceCells()
+                nbrFaceCells
+                //cpp.neighbPatch().faceCells()
             )
         )
     );
@@ -209,6 +217,7 @@ Foam::cyclicACMIFvPatchField<Type>::nonOverlapPatchField() const
             this->primitiveField()
         );
 
+    // WIP: Needs to re-direct nonOverlapPatchID to new patchId for assembly?
     return fld.boundaryField()[cyclicACMIPatch_.nonOverlapPatchID()];
 }
 
@@ -218,17 +227,24 @@ void Foam::cyclicACMIFvPatchField<Type>::updateInterfaceMatrix
 (
     solveScalarField& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const solveScalarField& psiInternal,
     const scalarField& coeffs,
     const direction cmpt,
     const Pstream::commsTypes
 ) const
 {
-    const cyclicACMIPolyPatch& cpp = cyclicACMIPatch_.cyclicACMIPatch();
+    // note: only applying coupled contribution
 
-    // Note: only applying coupled contribution
+//     const labelUList& nbrFaceCellsCoupled =
+//         lduAddr.patchAddr
+//         (
+//             cyclicACMIPatch_.cyclicACMIPatch().neighbPatchID()
+//         );
 
-    const labelUList& nbrFaceCellsCoupled = cpp.neighbPatch().faceCells();
+    const labelUList& nbrFaceCellsCoupled =
+        lduAddr.patchAddr(cyclicACMIPatch_.neighbPatchID());
 
     solveScalarField pnf(psiInternal, nbrFaceCellsCoupled);
 
@@ -237,7 +253,9 @@ void Foam::cyclicACMIFvPatchField<Type>::updateInterfaceMatrix
 
     pnf = cyclicACMIPatch_.interpolate(pnf);
 
-    this->addToInternalField(result, !add, coeffs, pnf);
+    const labelUList& faceCells = lduAddr.patchAddr(patchId);
+
+    this->addToInternalField(result, !add, faceCells, coeffs, pnf);
 }
 
 
@@ -246,16 +264,17 @@ void Foam::cyclicACMIFvPatchField<Type>::updateInterfaceMatrix
 (
     Field<Type>& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const Field<Type>& psiInternal,
     const scalarField& coeffs,
     const Pstream::commsTypes
 ) const
 {
-    const cyclicACMIPolyPatch& cpp = cyclicACMIPatch_.cyclicACMIPatch();
+    // note: only applying coupled contribution
 
-    // Note: only applying coupled contribution
-
-    const labelUList& nbrFaceCellsCoupled = cpp.neighbPatch().faceCells();
+    const labelUList& nbrFaceCellsCoupled =
+        lduAddr.patchAddr(cyclicACMIPatch_.neighbPatchID());
 
     Field<Type> pnf(psiInternal, nbrFaceCellsCoupled);
 
@@ -264,7 +283,9 @@ void Foam::cyclicACMIFvPatchField<Type>::updateInterfaceMatrix
 
     pnf = cyclicACMIPatch_.interpolate(pnf);
 
-    this->addToInternalField(result, !add, coeffs, pnf);
+    const labelUList& faceCells = lduAddr.patchAddr(patchId);
+
+    this->addToInternalField(result, !add, faceCells, coeffs, pnf);
 }
 
 
@@ -284,6 +305,141 @@ void Foam::cyclicACMIFvPatchField<Type>::manipulateMatrix
 }
 
 
+template<class Type>
+void Foam::cyclicACMIFvPatchField<Type>::manipulateMatrix
+(
+    fvMatrix<Type>& matrix,
+    const label mat,
+    const direction cmpt
+)
+{
+    if (this->cyclicACMIPatch().owner())
+    {
+        label index = this->patch().index();
+
+        const label globalPatchID =
+            matrix.lduMeshAssembly().patchLocalToGlobalMap()[mat][index];
+
+        const Field<scalar> intCoeffsCmpt
+        (
+            matrix.internalCoeffs()[globalPatchID].component(cmpt)
+        );
+
+        const Field<scalar> boundCoeffsCmpt
+        (
+            matrix.boundaryCoeffs()[globalPatchID].component(cmpt)
+        );
+
+        tmp<Field<scalar>> tintCoeffs(coeffs(matrix, intCoeffsCmpt, mat));
+        tmp<Field<scalar>> tbndCoeffs(coeffs(matrix, boundCoeffsCmpt, mat));
+        const Field<scalar>& intCoeffs = tintCoeffs.ref();
+        const Field<scalar>& bndCoeffs = tbndCoeffs.ref();
+
+        const labelUList& u = matrix.lduAddr().upperAddr();
+        const labelUList& l = matrix.lduAddr().lowerAddr();
+
+        label subFaceI = 0;
+
+        const labelList& faceMap =
+            matrix.lduMeshAssembly().faceBoundMap()[mat][index];
+
+        forAll (faceMap, j)
+        {
+            label globalFaceI = faceMap[j];
+
+            const scalar boundCorr = -bndCoeffs[subFaceI];
+            const scalar intCorr = -intCoeffs[subFaceI];
+
+            matrix.upper()[globalFaceI] += boundCorr;
+            matrix.diag()[u[globalFaceI]] -= intCorr;
+            matrix.diag()[l[globalFaceI]] -= boundCorr;
+
+            if (matrix.asymmetric())
+            {
+                matrix.lower()[globalFaceI] += intCorr;
+            }
+            subFaceI++;
+        }
+
+        // Set internalCoeffs and boundaryCoeffs in the assembly matrix
+        // on clyclicAMI patches to be used in the individual matrix by
+        // matrix.flux()
+        if (matrix.psi(mat).mesh().fluxRequired(this->internalField().name()))
+        {
+            matrix.internalCoeffs().set
+            (
+                globalPatchID, intCoeffs*pTraits<Type>::one
+            );
+            matrix.boundaryCoeffs().set
+            (
+                globalPatchID, bndCoeffs*pTraits<Type>::one
+            );
+
+            const label nbrPathID =
+                cyclicACMIPatch_.cyclicACMIPatch().neighbPatchID();
+
+            const label nbrGlobalPatchID =
+                matrix.lduMeshAssembly().patchLocalToGlobalMap()
+                [mat][nbrPathID];
+
+            matrix.internalCoeffs().set
+            (
+                nbrGlobalPatchID, intCoeffs*pTraits<Type>::one
+            );
+            matrix.boundaryCoeffs().set
+            (
+                nbrGlobalPatchID, bndCoeffs*pTraits<Type>::one
+            );
+        }
+    }
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Foam::scalar>>
+Foam::cyclicACMIFvPatchField<Type>::coeffs
+(
+    fvMatrix<Type>& matrix,
+    const Field<scalar>& coeffs,
+    const label mat
+) const
+{
+    const label index(this->patch().index());
+
+    const label nSubFaces
+    (
+        matrix.lduMeshAssembly().cellBoundMap()[mat][index].size()
+    );
+
+    Field<scalar> mapCoeffs(nSubFaces, Zero);
+
+    const scalarListList& srcWeight =
+        cyclicACMIPatch_.cyclicACMIPatch().AMI().srcWeights();
+
+    const scalarField& mask = cyclicACMIPatch_.cyclicACMIPatch().mask();
+
+    const scalar tol = cyclicACMIPolyPatch::tolerance();
+    label subFaceI = 0;
+    forAll(mask, faceI)
+    {
+        const scalarList& w = srcWeight[faceI];
+        for(label i=0; i<w.size(); i++)
+        {
+            if (mask[faceI] > tol)
+            {
+                const label localFaceId =
+                    matrix.lduMeshAssembly().facePatchFaceMap()
+                    [mat][index][subFaceI];
+                mapCoeffs[subFaceI] = w[i]*coeffs[localFaceId];
+            }
+            subFaceI++;
+        }
+    }
+
+    return tmp<Field<scalar>>(new Field<scalar>(mapCoeffs));
+}
+
+
 template<class Type>
 void Foam::cyclicACMIFvPatchField<Type>::updateCoeffs()
 {

src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.H

@@ -89,6 +89,15 @@ class cyclicACMIFvPatchField
             const Field<Type2>&
         ) const;
 
+         //- Return new matrix coeffs
+        tmp<Field<scalar>> coeffs
+        (
+            fvMatrix<Type>& matrix,
+            const Field<scalar>&,
+            const label
+        ) const;
+
+
 
 public:
 
@@ -206,6 +215,8 @@ public:
             (
                 solveScalarField& result,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const solveScalarField& psiInternal,
                 const scalarField& coeffs,
                 const direction cmpt,
@@ -217,6 +228,8 @@ public:
             (
                 Field<Type>&,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const Field<Type>&,
                 const scalarField&,
                 const Pstream::commsTypes commsType
@@ -225,6 +238,14 @@ public:
             //- Manipulate matrix
             virtual void manipulateMatrix(fvMatrix<Type>& matrix);
 
+            //- Manipulate matrix
+            virtual void manipulateMatrix
+            (
+                fvMatrix<Type>& m,
+                const label iMatrix,
+                const direction cmpt
+            );
+
             //- Update the coefficients associated with the patch field
             virtual void updateCoeffs();
 

src/finiteVolume/fields/fvPatchFields/constraint/cyclicAMI/cyclicAMIFvPatchField.C

@@ -26,6 +26,10 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+#include "fvMatrix.H"
+#include "volFields.H"
+//#include "cylicFvPatchField.H"
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class Type>
@@ -142,15 +146,20 @@ Foam::tmp<Foam::Field<Type>>
 Foam::cyclicAMIFvPatchField<Type>::patchNeighbourField() const
 {
     const Field<Type>& iField = this->primitiveField();
-    const labelUList& nbrFaceCells =
-        cyclicAMIPatch_.cyclicAMIPatch().neighbPatch().faceCells();
+
+    // By pass polyPatch to get nbrId. Instead use cyclicAMIFvPatch virtual
+    // neighbPatch()
+    const cyclicAMIFvPatch& neighbPatch = cyclicAMIPatch_.neighbPatch();
+    const labelUList& nbrFaceCells = neighbPatch.faceCells();
 
     Field<Type> pnf(iField, nbrFaceCells);
 
     tmp<Field<Type>> tpnf;
     if (cyclicAMIPatch_.applyLowWeightCorrection())
     {
-        tpnf = cyclicAMIPatch_.interpolate(pnf, this->patchInternalField()());
+        Field<Type> pnfInternal(iField, cyclicAMIPatch_.faceCells());
+
+        tpnf = cyclicAMIPatch_.interpolate(pnf, pnfInternal);
     }
     else
     {
@@ -188,6 +197,8 @@ void Foam::cyclicAMIFvPatchField<Type>::updateInterfaceMatrix
 (
     solveScalarField& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const solveScalarField& psiInternal,
     const scalarField& coeffs,
     const direction cmpt,
@@ -195,7 +206,7 @@ void Foam::cyclicAMIFvPatchField<Type>::updateInterfaceMatrix
 ) const
 {
     const labelUList& nbrFaceCells =
-        cyclicAMIPatch_.cyclicAMIPatch().neighbPatch().faceCells();
+        lduAddr.patchAddr(cyclicAMIPatch_.neighbPatchID());
 
     solveScalarField pnf(psiInternal, nbrFaceCells);
 
@@ -212,8 +223,10 @@ void Foam::cyclicAMIFvPatchField<Type>::updateInterfaceMatrix
         pnf = cyclicAMIPatch_.interpolate(pnf);
     }
 
+    const labelUList& faceCells = lduAddr.patchAddr(patchId);
+
     // Multiply the field by coefficients and add into the result
-    this->addToInternalField(result, !add, coeffs, pnf);
+    this->addToInternalField(result, !add, faceCells, coeffs, pnf);
 }
 
 
@@ -222,13 +235,15 @@ void Foam::cyclicAMIFvPatchField<Type>::updateInterfaceMatrix
 (
     Field<Type>& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const Field<Type>& psiInternal,
     const scalarField& coeffs,
     const Pstream::commsTypes
 ) const
 {
     const labelUList& nbrFaceCells =
-        cyclicAMIPatch_.cyclicAMIPatch().neighbPatch().faceCells();
+        lduAddr.patchAddr(cyclicAMIPatch_.neighbPatchID());
 
     Field<Type> pnf(psiInternal, nbrFaceCells);
 
@@ -245,8 +260,178 @@ void Foam::cyclicAMIFvPatchField<Type>::updateInterfaceMatrix
         pnf = cyclicAMIPatch_.interpolate(pnf);
     }
 
+    const labelUList& faceCells = lduAddr.patchAddr(patchId);
+
     // Multiply the field by coefficients and add into the result
-    this->addToInternalField(result, !add, coeffs, pnf);
+    this->addToInternalField(result, !add, faceCells, coeffs, pnf);
+}
+
+
+template<class Type>
+void Foam::cyclicAMIFvPatchField<Type>::manipulateMatrix
+(
+    fvMatrix<Type>& matrix,
+    const label mat,
+    const direction cmpt
+)
+{
+
+    if (this->cyclicAMIPatch().owner())
+    {
+        label index = this->patch().index();
+
+        const label globalPatchID =
+            matrix.lduMeshAssembly().patchLocalToGlobalMap()[mat][index];
+
+        const Field<scalar> intCoeffsCmpt
+        (
+            matrix.internalCoeffs()[globalPatchID].component(cmpt)
+        );
+
+        const Field<scalar> boundCoeffsCmpt
+        (
+            matrix.boundaryCoeffs()[globalPatchID].component(cmpt)
+        );
+
+        tmp<Field<scalar>> tintCoeffs(coeffs(matrix, intCoeffsCmpt, mat));
+        tmp<Field<scalar>> tbndCoeffs(coeffs(matrix, boundCoeffsCmpt, mat));
+        const Field<scalar>& intCoeffs = tintCoeffs.ref();
+        const Field<scalar>& bndCoeffs = tbndCoeffs.ref();
+
+        const labelUList& u = matrix.lduAddr().upperAddr();
+        const labelUList& l = matrix.lduAddr().lowerAddr();
+
+        label subFaceI = 0;
+
+        const labelList& faceMap =
+            matrix.lduMeshAssembly().faceBoundMap()[mat][index];
+
+        forAll (faceMap, j)
+        {
+            label globalFaceI = faceMap[j];
+
+            const scalar boundCorr = -bndCoeffs[subFaceI];
+            const scalar intCorr = -intCoeffs[subFaceI];
+
+            matrix.upper()[globalFaceI] += boundCorr;
+            matrix.diag()[u[globalFaceI]] -= intCorr;
+            matrix.diag()[l[globalFaceI]] -= boundCorr;
+
+            if (matrix.asymmetric())
+            {
+                matrix.lower()[globalFaceI] += intCorr;
+            }
+            subFaceI++;
+        }
+
+        // Set internalCoeffs and boundaryCoeffs in the assembly matrix
+        // on clyclicAMI patches to be used in the individual matrix by
+        // matrix.flux()
+        if (matrix.psi(mat).mesh().fluxRequired(this->internalField().name()))
+        {
+            matrix.internalCoeffs().set
+            (
+                globalPatchID, intCoeffs*pTraits<Type>::one
+            );
+            matrix.boundaryCoeffs().set
+            (
+                globalPatchID, bndCoeffs*pTraits<Type>::one
+            );
+
+            const label nbrPathID =
+                cyclicAMIPatch_.cyclicAMIPatch().neighbPatchID();
+
+            const label nbrGlobalPatchID =
+                matrix.lduMeshAssembly().patchLocalToGlobalMap()
+                    [mat][nbrPathID];
+
+            matrix.internalCoeffs().set
+            (
+                nbrGlobalPatchID, intCoeffs*pTraits<Type>::one
+            );
+            matrix.boundaryCoeffs().set
+            (
+                nbrGlobalPatchID, bndCoeffs*pTraits<Type>::one
+            );
+        }
+    }
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Foam::scalar>>
+Foam::cyclicAMIFvPatchField<Type>::coeffs
+(
+    fvMatrix<Type>& matrix,
+    const Field<scalar>& coeffs,
+    const label mat
+) const
+{
+    const label index(this->patch().index());
+
+    const label nSubFaces
+    (
+        matrix.lduMeshAssembly().cellBoundMap()[mat][index].size()
+    );
+
+    Field<scalar> mapCoeffs(nSubFaces, Zero);
+
+    const scalarListList& srcWeight =
+        cyclicAMIPatch_.cyclicAMIPatch().AMI().srcWeights();
+
+    label subFaceI = 0;
+    forAll(*this, faceI)
+    {
+        const scalarList& w = srcWeight[faceI];
+        for(label i=0; i<w.size(); i++)
+        {
+            const label localFaceId =
+                matrix.lduMeshAssembly().facePatchFaceMap()[mat][index][subFaceI];
+            mapCoeffs[subFaceI] = w[i]*coeffs[localFaceId];
+            subFaceI++;
+        }
+    }
+
+    return tmp<Field<scalar>>(new Field<scalar>(mapCoeffs));
+}
+
+
+template<class Type>
+template<class Type2>
+void Foam::cyclicAMIFvPatchField<Type>::collectStencilData
+(
+    const refPtr<mapDistribute>& mapPtr,
+    const labelListList& stencil,
+    const Type2& data,
+    List<Type2>& expandedData
+)
+{
+    expandedData.setSize(stencil.size());
+    if (mapPtr.valid())
+    {
+        Type2 work(data);
+        mapPtr().distribute(work);
+
+        forAll(stencil, facei)
+        {
+            const labelList& slots = stencil[facei];
+            expandedData[facei].append
+            (
+                UIndirectList<typename Type2::value_type>(work, slots)
+            );
+        }
+    }
+    else
+    {
+        forAll(stencil, facei)
+        {
+            const labelList& slots = stencil[facei];
+            expandedData[facei].append
+            (
+                UIndirectList<typename Type2::value_type>(data, slots)
+            );
+        }
+    }
 }
 
 
@@ -259,3 +444,4 @@ void Foam::cyclicAMIFvPatchField<Type>::write(Ostream& os) const
 
 
 // ************************************************************************* //
+

src/finiteVolume/fields/fvPatchFields/constraint/cyclicAMI/cyclicAMIFvPatchField.H

@@ -60,6 +60,7 @@ SourceFiles
 #ifndef cyclicAMIFvPatchField_H
 #define cyclicAMIFvPatchField_H
 
+
 #include "coupledFvPatchField.H"
 #include "cyclicAMILduInterfaceField.H"
 #include "cyclicAMIFvPatch.H"
@@ -95,6 +96,25 @@ class cyclicAMIFvPatchField
         ) const;
 
 
+        //- Return new matrix coeffs
+        tmp<Field<scalar>> coeffs
+        (
+            fvMatrix<Type>& matrix,
+            const Field<scalar>&,
+            const label
+        ) const;
+
+
+        template<class Type2>
+        void collectStencilData
+        (
+            const refPtr<mapDistribute>& mapPtr,
+            const labelListList& stencil,
+            const Type2& data,
+            List<Type2>& expandedData
+        );
+
+
 public:
 
     //- Runtime type information
@@ -187,6 +207,8 @@ public:
             (
                 solveScalarField& result,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const solveScalarField& psiInternal,
                 const scalarField& coeffs,
                 const direction cmpt,
@@ -198,12 +220,23 @@ public:
             (
                 Field<Type>&,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const Field<Type>&,
                 const scalarField&,
                 const Pstream::commsTypes commsType
             ) const;
 
 
+            //- Manipulate matrix
+            virtual void manipulateMatrix
+            (
+                fvMatrix<Type>& m,
+                const label iMatrix,
+                const direction cmpt
+            );
+
+
         // Cyclic AMI coupled interface functions
 
             //- Does the patch field perform the transformation

src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.C

@@ -101,8 +101,8 @@ Foam::jumpCyclicFvPatchField<Type>::patchNeighbourField() const
     const labelUList& nbrFaceCells =
         this->cyclicPatch().neighbFvPatch().faceCells();
 
-    tmp<Field<Type>> tpnf(new Field<Type>(this->size()));
-    Field<Type>& pnf = tpnf.ref();
+    auto tpnf = tmp<Field<Type>>::New(this->size());
+    auto& pnf = tpnf.ref();
 
     Field<Type> jf(this->jump());
     if (!this->cyclicPatch().owner())
@@ -137,6 +137,8 @@ void Foam::jumpCyclicFvPatchField<Type>::updateInterfaceMatrix
 (
     solveScalarField& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const solveScalarField& psiInternal,
     const scalarField& coeffs,
     const direction cmpt,
@@ -152,6 +154,8 @@ void Foam::jumpCyclicFvPatchField<Type>::updateInterfaceMatrix
 (
     Field<Type>& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const Field<Type>& psiInternal,
     const scalarField& coeffs,
     const Pstream::commsTypes
@@ -160,7 +164,10 @@ void Foam::jumpCyclicFvPatchField<Type>::updateInterfaceMatrix
     Field<Type> pnf(this->size());
 
     const labelUList& nbrFaceCells =
-        this->cyclicPatch().neighbFvPatch().faceCells();
+        lduAddr.patchAddr
+        (
+            this->cyclicPatch().neighbPatchID()
+        );
 
     // only apply jump to original field
     if (&psiInternal == &this->primitiveField())
@@ -188,8 +195,10 @@ void Foam::jumpCyclicFvPatchField<Type>::updateInterfaceMatrix
     // Transform according to the transformation tensors
     this->transformCoupleField(pnf);
 
+    const labelUList& faceCells = lduAddr.patchAddr(patchId);
+
     // Multiply the field by coefficients and add into the result
-    this->addToInternalField(result, !add, coeffs, pnf);
+    this->addToInternalField(result, !add, faceCells, coeffs, pnf);
 }
 
 

src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.H

@@ -133,6 +133,8 @@ public:
             (
                 solveScalarField& result,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const solveScalarField& psiInternal,
                 const scalarField& coeffs,
                 const direction cmpt,
@@ -144,6 +146,8 @@ public:
             (
                 Field<Type>&,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const Field<Type>&,
                 const scalarField&,
                 const Pstream::commsTypes commsType
@@ -156,6 +160,8 @@ void jumpCyclicFvPatchField<scalar>::updateInterfaceMatrix
 (
     solveScalarField& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const solveScalarField& psiInternal,
     const scalarField& coeffs,
     const direction cmpt,
@@ -168,6 +174,8 @@ void jumpCyclicFvPatchField<vector>::updateInterfaceMatrix
 (
     solveScalarField& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const solveScalarField& psiInternal,
     const scalarField& coeffs,
     const direction cmpt,

src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchFields.C

@@ -44,6 +44,8 @@ void Foam::jumpCyclicFvPatchField<Foam::scalar>::updateInterfaceMatrix
 (
     solveScalarField& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const solveScalarField& psiInternal,
     const scalarField& coeffs,
     const direction cmpt,
@@ -53,7 +55,10 @@ void Foam::jumpCyclicFvPatchField<Foam::scalar>::updateInterfaceMatrix
     solveScalarField pnf(this->size());
 
     const labelUList& nbrFaceCells =
-        this->cyclicPatch().neighbFvPatch().faceCells();
+        lduAddr.patchAddr
+        (
+            this->cyclicPatch().neighbPatchID()
+        );
 
     // only apply jump to original field
     if
@@ -85,8 +90,10 @@ void Foam::jumpCyclicFvPatchField<Foam::scalar>::updateInterfaceMatrix
     // Transform according to the transformation tensors
     this->transformCoupleField(pnf, cmpt);
 
+    const labelUList& faceCells = lduAddr.patchAddr(patchId);
+
     // Multiply the field by coefficients and add into the result
-    this->addToInternalField(result, !add, coeffs, pnf);
+    this->addToInternalField(result, !add, faceCells, coeffs, pnf);
 }
 
 
@@ -95,6 +102,8 @@ void Foam::jumpCyclicFvPatchField<Foam::vector>::updateInterfaceMatrix
 (
     solveScalarField& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const solveScalarField& psiInternal,
     const scalarField& coeffs,
     const direction cmpt,
@@ -104,7 +113,10 @@ void Foam::jumpCyclicFvPatchField<Foam::vector>::updateInterfaceMatrix
     solveScalarField pnf(this->size());
 
     const labelUList& nbrFaceCells =
-        this->cyclicPatch().neighbFvPatch().faceCells();
+        lduAddr.patchAddr
+        (
+            this->cyclicPatch().neighbPatchID()
+        );
 
     const Field<vector>& iField = this->primitiveField();
 
@@ -140,8 +152,10 @@ void Foam::jumpCyclicFvPatchField<Foam::vector>::updateInterfaceMatrix
     // Transform according to the transformation tensors
     this->transformCoupleField(pnf, cmpt);
 
+    const labelUList& faceCells = lduAddr.patchAddr(patchId);
+
     // Multiply the field by coefficients and add into the result
-    this->addToInternalField(result, !add, coeffs, pnf);
+    this->addToInternalField(result, !add, faceCells, coeffs, pnf);
 }
 
 // ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchField.C

@@ -138,6 +138,8 @@ void Foam::jumpCyclicAMIFvPatchField<Type>::updateInterfaceMatrix
 (
     solveScalarField& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const solveScalarField& psiInternal,
     const scalarField& coeffs,
     const direction cmpt,
@@ -153,13 +155,18 @@ void Foam::jumpCyclicAMIFvPatchField<Type>::updateInterfaceMatrix
 (
     Field<Type>& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const Field<Type>& psiInternal,
     const scalarField& coeffs,
     const Pstream::commsTypes
 ) const
 {
     const labelUList& nbrFaceCells =
-        this->cyclicAMIPatch().cyclicAMIPatch().neighbPatch().faceCells();
+        lduAddr.patchAddr
+        (
+            this->cyclicAMIPatch().neighbPatchID()
+        );
 
     Field<Type> pnf(psiInternal, nbrFaceCells);
 
@@ -193,8 +200,10 @@ void Foam::jumpCyclicAMIFvPatchField<Type>::updateInterfaceMatrix
     // Transform according to the transformation tensors
     this->transformCoupleField(pnf);
 
+    const labelUList& faceCells = lduAddr.patchAddr(patchId);
+
     // Multiply the field by coefficients and add into the result
-    this->addToInternalField(result, !add, coeffs, pnf);
+    this->addToInternalField(result, !add,faceCells, coeffs, pnf);
 }
 
 

src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchField.H

@@ -135,6 +135,8 @@ public:
             (
                 solveScalarField& result,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const solveScalarField& psiInternal,
                 const scalarField& coeffs,
                 const direction cmpt,
@@ -146,6 +148,8 @@ public:
             (
                 Field<Type>&,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const Field<Type>&,
                 const scalarField&,
                 const Pstream::commsTypes commsType
@@ -159,6 +163,8 @@ void jumpCyclicAMIFvPatchField<scalar>::updateInterfaceMatrix
 (
     solveScalarField& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const solveScalarField& psiInternal,
     const scalarField& coeffs,
     const direction cmpt,

src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchFields.C

@@ -46,6 +46,8 @@ void Foam::jumpCyclicAMIFvPatchField<scalar>::updateInterfaceMatrix
 (
     solveScalarField& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const solveScalarField& psiInternal,
     const scalarField& coeffs,
     const direction cmpt,
@@ -53,7 +55,10 @@ void Foam::jumpCyclicAMIFvPatchField<scalar>::updateInterfaceMatrix
 ) const
 {
     const labelUList& nbrFaceCells =
-        this->cyclicAMIPatch().cyclicAMIPatch().neighbPatch().faceCells();
+        lduAddr.patchAddr
+        (
+            this->cyclicAMIPatch().neighbPatchID()
+        );
 
     solveScalarField pnf(psiInternal, nbrFaceCells);
 
@@ -83,8 +88,10 @@ void Foam::jumpCyclicAMIFvPatchField<scalar>::updateInterfaceMatrix
     // Transform according to the transformation tensors
     this->transformCoupleField(pnf, cmpt);
 
+    const labelUList& faceCells = lduAddr.patchAddr(patchId);
+
     // Multiply the field by coefficients and add into the result
-    this->addToInternalField(result, !add, coeffs, pnf);
+    this->addToInternalField(result, !add, faceCells, coeffs, pnf);
 }
 
 

src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchField.C

@@ -313,13 +313,23 @@ void Foam::processorFvPatchField<Type>::initInterfaceMatrixUpdate
 (
     solveScalarField&,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const solveScalarField& psiInternal,
     const scalarField&,
     const direction,
     const Pstream::commsTypes commsType
 ) const
 {
-    this->patch().patchInternalField(psiInternal, scalarSendBuf_);
+    //this->patch().patchInternalField(psiInternal, scalarSendBuf_);
+
+    const labelUList& faceCells = lduAddr.patchAddr(patchId);
+
+    scalarSendBuf_.setSize(this->patch().size());
+    forAll(scalarSendBuf_, facei)
+    {
+        scalarSendBuf_[facei] = psiInternal[faceCells[facei]];
+    }
 
     if
     (
@@ -374,6 +384,8 @@ void Foam::processorFvPatchField<Type>::updateInterfaceMatrix
 (
     solveScalarField& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const solveScalarField&,
     const scalarField& coeffs,
     const direction cmpt,
@@ -385,6 +397,8 @@ void Foam::processorFvPatchField<Type>::updateInterfaceMatrix
         return;
     }
 
+    const labelUList& faceCells = lduAddr.patchAddr(patchId);
+
     if
     (
         commsType == Pstream::commsTypes::nonBlocking
@@ -412,7 +426,14 @@ void Foam::processorFvPatchField<Type>::updateInterfaceMatrix
         }
 
         // Multiply the field by coefficients and add into the result
-        this->addToInternalField(result, !add, coeffs, scalarReceiveBuf_);
+        this->addToInternalField
+        (
+            result,
+            !add,
+            faceCells,
+            coeffs,
+            scalarReceiveBuf_
+        );
     }
     else
     {
@@ -432,7 +453,7 @@ void Foam::processorFvPatchField<Type>::updateInterfaceMatrix
         }
 
         // Multiply the field by coefficients and add into the result
-        this->addToInternalField(result, !add, coeffs, pnf);
+        this->addToInternalField(result, !add, faceCells, coeffs, pnf);
     }
 
     const_cast<processorFvPatchField<Type>&>(*this).updatedMatrix() = true;
@@ -444,12 +465,21 @@ void Foam::processorFvPatchField<Type>::initInterfaceMatrixUpdate
 (
     Field<Type>&,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const Field<Type>& psiInternal,
     const scalarField&,
     const Pstream::commsTypes commsType
 ) const
 {
-    this->patch().patchInternalField(psiInternal, sendBuf_);
+    sendBuf_.setSize(this->patch().size());
+
+    const labelUList& faceCells = lduAddr.patchAddr(patchId);
+
+    forAll(sendBuf_, facei)
+    {
+        sendBuf_[facei] = psiInternal[faceCells[facei]];
+    }
 
     if
     (
@@ -504,6 +534,8 @@ void Foam::processorFvPatchField<Type>::updateInterfaceMatrix
 (
     Field<Type>& result,
     const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
     const Field<Type>&,
     const scalarField& coeffs,
     const Pstream::commsTypes commsType
@@ -514,6 +546,8 @@ void Foam::processorFvPatchField<Type>::updateInterfaceMatrix
         return;
     }
 
+    const labelUList& faceCells = lduAddr.patchAddr(patchId);
+
     if
     (
         commsType == Pstream::commsTypes::nonBlocking
@@ -539,7 +573,7 @@ void Foam::processorFvPatchField<Type>::updateInterfaceMatrix
         transformCoupleField(receiveBuf_);
 
         // Multiply the field by coefficients and add into the result
-        this->addToInternalField(result, !add, coeffs, receiveBuf_);
+        this->addToInternalField(result, !add, faceCells, coeffs, receiveBuf_);
     }
     else
     {
@@ -552,7 +586,7 @@ void Foam::processorFvPatchField<Type>::updateInterfaceMatrix
         transformCoupleField(pnf);
 
         // Multiply the field by coefficients and add into the result
-        this->addToInternalField(result, !add, coeffs, pnf);
+        this->addToInternalField(result, !add, faceCells, coeffs, pnf);
     }
 
     const_cast<processorFvPatchField<Type>&>(*this).updatedMatrix() = true;

src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchField.H

@@ -214,6 +214,8 @@ public:
             (
                 solveScalarField& result,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const solveScalarField& psiInternal,
                 const scalarField& coeffs,
                 const direction cmpt,
@@ -225,6 +227,8 @@ public:
             (
                 solveScalarField& result,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const solveScalarField& psiInternal,
                 const scalarField& coeffs,
                 const direction cmpt,
@@ -236,6 +240,8 @@ public:
             (
                 Field<Type>& result,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const Field<Type>& psiInternal,
                 const scalarField& coeffs,
                 const Pstream::commsTypes commsType
@@ -246,10 +252,13 @@ public:
             (
                 Field<Type>& result,
                 const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
                 const Field<Type>& psiInternal,
                 const scalarField& coeffs,
                 const Pstream::commsTypes commsType
             ) const;
+            
 
 
         //- Processor coupled interface functions

src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedMixedFieldFvPatchField/mappedMixedFieldFvPatchFields.C

@@ -2,8 +2,10 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
+    \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.

src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedMixedFieldFvPatchField/mappedMixedFieldFvPatchFields.H

@@ -2,8 +2,10 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
+    \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.

src/finiteVolume/fields/fvPatchFields/derived/mappedMixed/mappedMixedFvPatchFields.C

@@ -2,8 +2,10 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2020 OpenCFD Ltd.
+    \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.

src/finiteVolume/fields/fvPatchFields/derived/mappedMixed/mappedMixedFvPatchFields.H

@@ -2,8 +2,10 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
+    \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.

src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.C

@@ -46,7 +46,8 @@ Foam::fvPatchField<Type>::fvPatchField
     internalField_(iF),
     updated_(false),
     manipulatedMatrix_(false),
-    patchType_(word::null)
+    patchType_(word::null),
+    useImplicit_(false)
 {}
 
 
@@ -63,7 +64,8 @@ Foam::fvPatchField<Type>::fvPatchField
     internalField_(iF),
     updated_(false),
     manipulatedMatrix_(false),
-    patchType_(word::null)
+    patchType_(word::null),
+    useImplicit_(false)
 {}
 
 
@@ -80,7 +82,8 @@ Foam::fvPatchField<Type>::fvPatchField
     internalField_(iF),
     updated_(false),
     manipulatedMatrix_(false),
-    patchType_(patchType)
+    patchType_(patchType),
+    useImplicit_(false)
 {}
 
 
@@ -97,7 +100,8 @@ Foam::fvPatchField<Type>::fvPatchField
     internalField_(iF),
     updated_(false),
     manipulatedMatrix_(false),
-    patchType_(word::null)
+    patchType_(word::null),
+    useImplicit_(false)
 {}
 
 
@@ -115,7 +119,8 @@ Foam::fvPatchField<Type>::fvPatchField
     internalField_(iF),
     updated_(false),
     manipulatedMatrix_(false),
-    patchType_(dict.getOrDefault<word>("patchType", word::null))
+    patchType_(dict.getOrDefault<word>("patchType", word::null)),
+    useImplicit_(dict.getOrDefault<bool>("useImplicit", false))
 {
     if (valueRequired)
     {
@@ -150,7 +155,8 @@ Foam::fvPatchField<Type>::fvPatchField
     internalField_(iF),
     updated_(false),
     manipulatedMatrix_(false),
-    patchType_(ptf.patchType_)
+    patchType_(ptf.patchType_),
+    useImplicit_(ptf.useImplicit_)
 {
     // For unmapped faces set to internal field value (zero-gradient)
     if (notNull(iF) && mapper.hasUnmapped())
@@ -172,7 +178,8 @@ Foam::fvPatchField<Type>::fvPatchField
     internalField_(ptf.internalField_),
     updated_(false),
     manipulatedMatrix_(false),
-    patchType_(ptf.patchType_)
+    patchType_(ptf.patchType_),
+    useImplicit_(ptf.useImplicit_)
 {}
 
 
@@ -188,7 +195,8 @@ Foam::fvPatchField<Type>::fvPatchField
     internalField_(iF),
     updated_(false),
     manipulatedMatrix_(false),
-    patchType_(ptf.patchType_)
+    patchType_(ptf.patchType_),
+    useImplicit_(ptf.useImplicit_)
 {}
 
 
@@ -360,10 +368,23 @@ void Foam::fvPatchField<Type>::manipulateMatrix
 }
 
 
+template<class Type>
+void Foam::fvPatchField<Type>::manipulateMatrix
+(
+    fvMatrix<Type>& matrix,
+    const label iMatrix,
+    const direction cmp
+)
+{
+    manipulatedMatrix_ = true;
+}
+
+
 template<class Type>
 void Foam::fvPatchField<Type>::write(Ostream& os) const
 {
     os.writeEntry("type", type());
+    os.writeEntryIfDifferent<bool>("useImplicit", false, useImplicit_);
 
     if (patchType_.size())
     {
@@ -585,6 +606,6 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const fvPatchField<Type>& ptf)
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-    #include "fvPatchFieldNew.C"
+#include "fvPatchFieldNew.C"
 
 // ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H

@@ -57,24 +57,17 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of classes
+// Forward Declarations
 
 class objectRegistry;
 class dictionary;
 class fvPatchFieldMapper;
 class volMesh;
 
+template<class Type> class fvPatchField;
+template<class Type> class calculatedFvPatchField;
+template<class Type> class fvMatrix;
 
-// Forward declaration of friend functions and operators
-
-template<class Type>
-class fvPatchField;
-
-template<class Type>
-class calculatedFvPatchField;
-
-template<class Type>
-class fvMatrix;
 
 template<class Type>
 Ostream& operator<<(Ostream&, const fvPatchField<Type>&);
@@ -110,6 +103,9 @@ class fvPatchField
         //  patch type as 'patchType'
         word patchType_;
 
+        //- Use implicit formulation
+        bool useImplicit_;
+
 
 public:
 
@@ -306,7 +302,22 @@ public:
     virtual ~fvPatchField<Type>() = default;
 
 
-    // Member functions
+    // Member Functions
+
+        // Implicit Functions
+
+            //- Use implicit formulation for coupled patches only
+            bool useImplicit() const
+            {
+                return useImplicit_;
+            }
+
+            //- Non-const access to useImplicit
+            bool& useImplicit()
+            {
+                return useImplicit_;
+            }
+
 
         // Attributes
 
@@ -382,7 +393,6 @@ public:
                 return manipulatedMatrix_;
             }
 
-
         // Mapping functions
 
             //- Map (and resize as needed) from self given a mapping object
@@ -527,6 +537,14 @@ public:
                 const scalarField& weights
             );
 
+            //- Manipulate fvMatrix
+            virtual void manipulateMatrix
+            (
+                fvMatrix<Type>& matrix,
+                const label iMatrix,
+                const direction cmp
+            );
+
 
         // I-O
 

src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -35,6 +35,12 @@ License
 #include "UniformList.H"
 #include "demandDrivenData.H"
 
+#include "cyclicFvPatchField.H"
+#include "cyclicAMIFvPatchField.H"
+#include "cyclicACMIFvPatchField.H"
+
+#include "processorLduInterfaceField.H"
+
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 template<class Type>
@@ -120,14 +126,24 @@ void Foam::fvMatrix<Type>::addBoundaryDiag
     const direction solveCmpt
 ) const
 {
-    forAll(internalCoeffs_, patchi)
+    for (label fieldi = 0; fieldi < nMatrices(); fieldi++)
     {
-        addToInternalField
-        (
-            lduAddr().patchAddr(patchi),
-            internalCoeffs_[patchi].component(solveCmpt),
-            diag
-        );
+        const auto& bpsi = this->psi(fieldi).boundaryField();
+
+        forAll(bpsi, ptfi)
+        {
+            const label patchi = globalPatchID(fieldi, ptfi);
+
+            if (patchi != -1)
+            {
+                addToInternalField
+                (
+                    lduAddr().patchAddr(patchi),
+                    internalCoeffs_[patchi].component(solveCmpt),
+                    diag
+                );
+            }
+        }
     }
 }
 
@@ -135,14 +151,23 @@ void Foam::fvMatrix<Type>::addBoundaryDiag
 template<class Type>
 void Foam::fvMatrix<Type>::addCmptAvBoundaryDiag(scalarField& diag) const
 {
-    forAll(internalCoeffs_, patchi)
+    for (label fieldi = 0; fieldi < nMatrices(); fieldi++)
     {
-        addToInternalField
-        (
-            lduAddr().patchAddr(patchi),
-            cmptAv(internalCoeffs_[patchi]),
-            diag
-        );
+        const auto& bpsi = this->psi(fieldi).boundaryField();
+
+        forAll(bpsi, ptfi)
+        {
+            const label patchi = globalPatchID(fieldi, ptfi);
+            if (patchi != -1)
+            {
+                addToInternalField
+                (
+                    lduAddr().patchAddr(patchi),
+                    cmptAv(internalCoeffs_[patchi]),
+                    diag
+                );
+            }
+        }
     }
 }
 
@@ -154,25 +179,42 @@ void Foam::fvMatrix<Type>::addBoundarySource
     const bool couples
 ) const
 {
-    forAll(psi_.boundaryField(), patchi)
+    for (label fieldi = 0; fieldi < nMatrices(); fieldi++)
     {
-        const fvPatchField<Type>& ptf = psi_.boundaryField()[patchi];
-        const Field<Type>& pbc = boundaryCoeffs_[patchi];
+        const auto& bpsi = this->psi(fieldi).boundaryField();
 
-        if (!ptf.coupled())
+        forAll(bpsi, ptfi)
         {
-            addToInternalField(lduAddr().patchAddr(patchi), pbc, source);
-        }
-        else if (couples)
-        {
-            const tmp<Field<Type>> tpnf = ptf.patchNeighbourField();
-            const Field<Type>& pnf = tpnf();
+            const fvPatchField<Type>& ptf = bpsi[ptfi];
 
-            const labelUList& addr = lduAddr().patchAddr(patchi);
+            const label patchi = globalPatchID(fieldi, ptfi);
 
-            forAll(addr, facei)
+            if (patchi != -1)
             {
-                source[addr[facei]] += cmptMultiply(pbc[facei], pnf[facei]);
+                const Field<Type>& pbc = boundaryCoeffs_[patchi];
+
+                if (!ptf.coupled())
+                {
+                    addToInternalField
+                    (
+                        lduAddr().patchAddr(patchi),
+                        pbc,
+                        source
+                    );
+                }
+                else if (couples)
+                {
+                    const tmp<Field<Type>> tpnf = ptf.patchNeighbourField();
+                    const Field<Type>& pnf = tpnf();
+
+                    const labelUList& addr = lduAddr().patchAddr(patchi);
+
+                    forAll(addr, facei)
+                    {
+                        source[addr[facei]] +=
+                            cmptMultiply(pbc[facei], pnf[facei]);
+                    }
+                }
             }
         }
     }
@@ -261,8 +303,31 @@ void Foam::fvMatrix<Type>::setValuesFromList
 }
 
 
+template<class Type>
+Foam::label Foam::fvMatrix<Type>::checkImplicit(const label fieldI)
+{
+    const auto& bpsi = this->psi(fieldI).boundaryField();
+
+    label id = -1;
+    forAll (bpsi, patchI)
+    {
+        if (bpsi[patchI].useImplicit())
+        {
+            id += (label(2) << patchI);
+        }
+    }
+    if (id > 0)
+    {
+        lduAssemblyName_ = word("lduAssembly") + name(id);
+        useImplicit_ = true;
+    }
+    return id;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+
 template<class Type>
 Foam::fvMatrix<Type>::fvMatrix
 (
@@ -272,6 +337,9 @@ Foam::fvMatrix<Type>::fvMatrix
 :
     lduMatrix(psi.mesh()),
     psi_(psi),
+    useImplicit_(false),
+    lduAssemblyName_(word::null),
+    nMatrix_(0),
     dimensions_(ds),
     source_(psi.size(), Zero),
     internalCoeffs_(psi.mesh().boundary().size()),
@@ -281,7 +349,8 @@ Foam::fvMatrix<Type>::fvMatrix
     DebugInFunction
         << "Constructing fvMatrix<Type> for field " << psi_.name() << endl;
 
-    // Initialise coupling coefficients
+    checkImplicit();
+
     forAll(psi.mesh().boundary(), patchi)
     {
         internalCoeffs_.set
@@ -305,11 +374,8 @@ Foam::fvMatrix<Type>::fvMatrix
         );
     }
 
-    // Update the boundary coefficients of psi without changing its event No.
-    auto& psiRef =
-        const_cast<GeometricField<Type, fvPatchField, volMesh>&>(psi_);
-
-    const label currentStatePsi = psiRef.eventNo();
+    auto& psiRef = this->psi(0);
+    label currentStatePsi = psiRef.eventNo();
     psiRef.boundaryFieldRef().updateCoeffs();
     psiRef.eventNo() = currentStatePsi;
 }
@@ -321,6 +387,9 @@ Foam::fvMatrix<Type>::fvMatrix(const fvMatrix<Type>& fvm)
     refCount(),
     lduMatrix(fvm),
     psi_(fvm.psi_),
+    useImplicit_(fvm.useImplicit_),
+    lduAssemblyName_(fvm.lduAssemblyName_),
+    nMatrix_(fvm.nMatrix_),
     dimensions_(fvm.dimensions_),
     source_(fvm.source_),
     internalCoeffs_(fvm.internalCoeffs_),
@@ -350,6 +419,9 @@ Foam::fvMatrix<Type>::fvMatrix(const tmp<fvMatrix<Type>>& tfvm)
         tfvm.isTmp()
     ),
     psi_(tfvm().psi_),
+    useImplicit_(tfvm().useImplicit_),
+    lduAssemblyName_(tfvm().lduAssemblyName_),
+    nMatrix_(tfvm().nMatrix_),
     dimensions_(tfvm().dimensions_),
     source_
     (
@@ -401,15 +473,21 @@ Foam::fvMatrix<Type>::fvMatrix
 :
     lduMatrix(psi.mesh()),
     psi_(psi),
+    useImplicit_(false),
+    lduAssemblyName_(word::null),
+    nMatrix_(0),
     dimensions_(is),
     source_(is),
     internalCoeffs_(psi.mesh().boundary().size()),
     boundaryCoeffs_(psi.mesh().boundary().size()),
     faceFluxCorrectionPtr_(nullptr)
 {
+
     DebugInFunction
         << "Constructing fvMatrix<Type> for field " << psi_.name() << endl;
 
+    checkImplicit();
+
     // Initialise coupling coefficients
     forAll(psi.mesh().boundary(), patchi)
     {
@@ -455,11 +533,519 @@ Foam::fvMatrix<Type>::~fvMatrix()
         << "Destroying fvMatrix<Type> for field " << psi_.name() << endl;
 
     deleteDemandDrivenData(faceFluxCorrectionPtr_);
+    subMatrices_.clear();
 }
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class Type>
+void Foam::fvMatrix<Type>::setInterfaces
+(
+    lduInterfaceFieldPtrsList& interfaces,
+    PtrDynList<lduInterfaceField>& newInterfaces
+)
+{
+    interfaces.setSize(internalCoeffs_.size());
+    for (label fieldi = 0; fieldi < nMatrices(); fieldi++)
+    {
+        const auto& bpsi = this->psi(fieldi).boundaryField();
+        lduInterfaceFieldPtrsList fieldInterfaces(bpsi.scalarInterfaces());
+
+        forAll (fieldInterfaces, patchi)
+        {
+            label globalPatchID = lduMeshPtr()->patchMap()[fieldi][patchi];
+
+            if (globalPatchID != -1)
+            {
+                if (fieldInterfaces.set(patchi))
+                {
+                    if (isA<cyclicLduInterfaceField>(bpsi[patchi]))
+                    {
+                        newInterfaces.append
+                        (
+                            new cyclicFvPatchField<Type>
+                            (
+                                refCast<const fvPatch>
+                                (
+                                    lduMeshPtr()->interfaces()[globalPatchID]
+                                ),
+                                bpsi[patchi].internalField()
+                            )
+                        );
+                        interfaces.set(globalPatchID, &newInterfaces.last());
+
+                    }
+                    else if (isA<cyclicAMILduInterfaceField>(bpsi[patchi]))
+                    {
+                        newInterfaces.append
+                        (
+                            new cyclicAMIFvPatchField<Type>
+                            (
+                                refCast<const fvPatch>
+                                (
+                                    lduMeshPtr()->interfaces()[globalPatchID]
+                                ),
+                                bpsi[patchi].internalField()
+                            )
+                        );
+                        interfaces.set(globalPatchID, &newInterfaces.last());
+                    }
+                    else if (isA<cyclicACMILduInterfaceField>(bpsi[patchi]))
+                    {
+                        newInterfaces.append
+                        (
+                            new cyclicACMIFvPatchField<Type>
+                            (
+                                refCast<const fvPatch>
+                                (
+                                    lduMeshPtr()->interfaces()[globalPatchID]
+                                ),
+                                bpsi[patchi].internalField()
+                            )
+                        );
+                        interfaces.set(globalPatchID, &newInterfaces.last());
+                    }
+                    else
+                    {
+                        interfaces.set(globalPatchID, &fieldInterfaces[patchi]);
+                    }
+                }
+            }
+        }
+    }
+}
+
+
+template<class Type>
+void Foam::fvMatrix<Type>::mapContributions
+(
+    label fieldi,
+    const FieldField<Field, Type>& fluxContrib,
+    FieldField<Field, Type>& contrib,
+    bool internal
+) const
+{
+    const lduPrimitiveMeshAssembly* ptr = lduMeshPtr();
+
+    const labelList& patchMap = ptr->patchMap()[fieldi];
+
+    forAll(contrib, patchi)
+    {
+        const label globalPtchId = patchMap[patchi];
+
+        if (globalPtchId != -1)
+        {
+            // Cache contrib before overwriting
+            const Field<Type> saveContrib(fluxContrib[globalPtchId]);
+            contrib[patchi].setSize(psi_.boundaryField()[patchi].size()),
+            contrib[patchi] = pTraits<Type>::zero;
+
+            if (internal)
+            {
+                contrib[patchi] =
+                    cmptMultiply
+                    (
+                        saveContrib,
+                        psi_.boundaryField()[patchi].patchInternalField()
+                    );
+            }
+            else
+            {
+                if (this->psi(fieldi).boundaryField()[patchi].coupled())
+                {
+                    contrib[patchi] =
+                        cmptMultiply
+                        (
+                            saveContrib,
+                            psi_.boundaryField()[patchi].patchNeighbourField()
+                        );
+                }
+            }
+        }
+        else if (globalPtchId == -1)
+        {
+            const polyPatch& pp =
+                this->psi(fieldi).mesh().boundaryMesh()[patchi];
+
+            if (pp.masterImplicit())
+            {
+                label virtualPatch =
+                    ptr->patchLocalToGlobalMap()[fieldi][patchi];
+
+                const label nbrPatchId = pp.neighbPolyPatchID();
+
+                // Copy contrib before overwriting
+                const Field<Type> saveContrib(fluxContrib[virtualPatch]);
+
+                Field<Type>& coeffs = contrib[patchi];
+                Field<Type>& nbrCoeffs = contrib[nbrPatchId];
+
+                coeffs.setSize(psi_.boundaryField()[patchi].size());
+                nbrCoeffs.setSize(psi_.boundaryField()[nbrPatchId].size());
+
+                coeffs = pTraits<Type>::zero;
+                nbrCoeffs = pTraits<Type>::zero;
+
+                // nrb cells
+                const labelList& nbrCellIds =
+                    ptr->cellBoundMap()[fieldi][patchi];
+
+                const labelList& cellIds =
+                    ptr->cellBoundMap()[fieldi][nbrPatchId];
+
+                const GeometricField<Type, fvPatchField, volMesh>& psi =
+                    this->psi(fieldi);
+
+                forAll(saveContrib, subFaceI)
+                {
+                    const label faceId =
+                        ptr->facePatchFaceMap()[fieldi][patchi][subFaceI];
+                    const label nbrFaceId =
+                        ptr->facePatchFaceMap()[fieldi][nbrPatchId][subFaceI];
+
+                    const label nbrCellId = nbrCellIds[subFaceI];
+                    const label cellId = cellIds[subFaceI];
+
+                    if (internal)
+                    {
+                        coeffs[faceId] +=
+                            cmptMultiply(saveContrib[subFaceI], psi[cellId]);
+
+                        nbrCoeffs[nbrFaceId] +=
+                            cmptMultiply(saveContrib[subFaceI], psi[nbrCellId]);
+                    }
+                    else //boundary
+                    {
+                        coeffs[faceId] +=
+                            cmptMultiply(saveContrib[subFaceI], psi[nbrCellId]);
+
+                        nbrCoeffs[nbrFaceId] +=
+                            cmptMultiply(saveContrib[subFaceI], psi[cellId]);
+                    }
+                }
+            }
+        }
+    }
+}
+
+
+template<class Type>
+void Foam::fvMatrix<Type>::setBounAndInterCoeffs()
+{
+    // If it is a multi-fvMatrix needs correct internalCoeffs and
+    // boundaryCoeffs size
+    if (nMatrix_ > 0)
+    {
+        label interfaceI(0);
+        for (label fieldi = 0; fieldi < nMatrices(); fieldi++)
+        {
+            const auto& psi = this->psi(fieldi);
+
+            forAll(psi.mesh().boundary(), patchi)
+            {
+                interfaceI++;
+            }
+        }
+        internalCoeffs_.setSize(interfaceI);
+        boundaryCoeffs_.setSize(interfaceI);
+
+        interfaceI = 0;
+        for (label fieldi = 0; fieldi < nMatrices(); fieldi++)
+        {
+            const auto& psi = this->psi(fieldi);
+
+            forAll(psi.mesh().boundary(), patchi)
+            {
+                internalCoeffs_.set
+                (
+                    interfaceI,
+                    new Field<Type>
+                    (
+                        psi.mesh().boundary()[patchi].size(),
+                        Zero
+                    )
+                );
+
+                boundaryCoeffs_.set
+                (
+                    interfaceI,
+                    new Field<Type>
+                    (
+                        psi.mesh().boundary()[patchi].size(),
+                        Zero
+                    )
+                );
+                interfaceI++;
+            }
+        }
+    }
+
+    for (label i=0; i < nMatrices(); ++i)
+    {
+        const auto& bpsi = this->psi(i).boundaryField();
+
+        // Cache to-be implicit internal/boundary
+        FieldField<Field, Type> boundary(bpsi.size());
+        FieldField<Field, Type> internal(bpsi.size());
+
+        label implicit = 0;
+        forAll(bpsi, patchI)
+        {
+            label globalPatchId = lduMeshPtr()->patchMap()[i][patchI];
+            if (globalPatchId == -1)
+            {
+                boundary.set
+                (
+                    implicit,
+                    matrix(i).boundaryCoeffs()[patchI].clone()
+                );
+                internal.set
+                (
+                    implicit,
+                    matrix(i).internalCoeffs()[patchI].clone()
+                );
+                implicit++;
+            }
+        }
+
+        // Update non-implicit patches (re-order)
+        forAll(bpsi, patchI)
+        {
+            label globalPatchId = lduMeshPtr()->patchMap()[i][patchI];
+            if (globalPatchId != -1)
+            {
+                if (matrix(i).internalCoeffs().set(patchI))
+                {
+                    internalCoeffs_.set
+                    (
+                        globalPatchId,
+                        matrix(i).internalCoeffs()[patchI].clone()
+                    );
+                }
+
+                if (matrix(i).boundaryCoeffs().set(patchI))
+                {
+                    boundaryCoeffs_.set
+                    (
+                        globalPatchId,
+                        matrix(i).boundaryCoeffs()[patchI].clone()
+                    );
+                }
+            }
+        }
+
+        // Store implicit patches at the end of the list
+        implicit = 0;
+        forAll(bpsi, patchI)
+        {
+            label globalPatchId = lduMeshPtr()->patchMap()[i][patchI];
+            if (globalPatchId == -1)
+            {
+                const label implicitPatchId =
+                     lduMeshPtr()->patchLocalToGlobalMap()[i][patchI];
+
+                internalCoeffs_.set
+                (
+                    implicitPatchId, internal[implicit].clone()
+                );
+                boundaryCoeffs_.set
+                (
+                    implicitPatchId, boundary[implicit].clone()
+                );
+
+                implicit++;
+            }
+        }
+    }
+
+//    forAll(internalCoeffs_, patchI)
+//    {
+//         DebugVar(patchI)
+//         DebugVar(internalCoeffs_[patchI])
+//         DebugVar(boundaryCoeffs_[patchI])
+//    }
+}
+
+
+template<class Type>
+void Foam::fvMatrix<Type>::manipulateMatrix(direction cmp)
+{
+    for (label i=0; i < nMatrices(); ++i)
+    {
+        forAll(psi(i).boundaryField(), patchI)
+        {
+            label globalPatchId = lduMeshPtr()->patchMap()[i][patchI];
+
+            if (globalPatchId == -1)
+            {
+                psi(i).boundaryFieldRef()[patchI].manipulateMatrix
+                (
+                    *this,
+                    i,
+                    cmp
+                );
+            }
+        }
+    }
+}
+
+
+template<class Type>
+void Foam::fvMatrix<Type>::transferFvMatrixCoeffs()
+{
+    const labelListList& faceMap = lduMeshPtr()->faceMap();
+    const labelList& cellMap = lduMeshPtr()->cellOffsets();
+
+    label newFaces = lduMeshPtr()->lduAddr().upperAddr().size();
+    label newCells = lduMeshPtr()->lduAddr().size();
+
+    scalarField lowerAssemb(newFaces, Zero);
+    scalarField upperAssemb(newFaces, Zero);
+    scalarField diagAssemb(newCells, Zero);
+    Field<Type> sourceAssemb(newCells, Zero);
+
+    bool asymmetricAssemby = false;
+    for (label i=0; i < nMatrices(); ++i)
+    {
+        if (matrix(i).asymmetric())
+        {
+            asymmetricAssemby = true;
+        }
+    }
+    // Move append contents into intermediate list
+    for (label i=0; i < nMatrices(); ++i)
+    {
+        if (asymmetricAssemby)
+        {
+            const scalarField lowerSub(matrix(i).lower());
+            forAll(lowerSub, facei)
+            {
+                lowerAssemb[faceMap[i][facei]] = lowerSub[facei];
+            }
+        }
+
+        const scalarField upperSub(matrix(i).upper());
+        const scalarField diagSub(matrix(i).diag());
+        const Field<Type> sourceSub(matrix(i).source());
+
+        forAll(upperSub, facei)
+        {
+            upperAssemb[faceMap[i][facei]] = upperSub[facei];
+        }
+
+        forAll(diagSub, celli)
+        {
+            const label globalCelli = cellMap[i] + celli;
+            diagAssemb[globalCelli] = diagSub[celli];
+            sourceAssemb[globalCelli] = sourceSub[celli];
+        }
+    }
+
+    if (asymmetricAssemby)
+    {
+        lower().setSize(newFaces, Zero);
+        lower() = lowerAssemb;
+    }
+    upper().setSize(newFaces, Zero);
+    upper() = upperAssemb;
+
+    diag().setSize(newCells, Zero);
+    diag() = diagAssemb;
+
+    source().setSize(newCells, Zero);
+    source() = sourceAssemb;
+}
+
+
+template<class Type>
+Foam::lduPrimitiveMeshAssembly* Foam::fvMatrix<Type>::lduMeshPtr()
+{
+    const lduPrimitiveMeshAssembly* lduAssemMeshPtr =
+        psi_.mesh().thisDb().objectRegistry::template findObject
+        <
+            lduPrimitiveMeshAssembly
+        > (lduAssemblyName_);
+
+    return const_cast<lduPrimitiveMeshAssembly*>(lduAssemMeshPtr);
+}
+
+
+template<class Type>
+const Foam::lduPrimitiveMeshAssembly* Foam::fvMatrix<Type>::lduMeshPtr() const
+{
+    return
+    (
+        psi_.mesh().thisDb().objectRegistry::template cfindObject
+        <
+            lduPrimitiveMeshAssembly
+        > (lduAssemblyName_)
+    );
+}
+
+
+template<class Type>
+void Foam::fvMatrix<Type>::createOrUpdateLduPrimitiveAssembly()
+{
+    lduPrimitiveMeshAssembly* ptr = lduMeshPtr();
+
+    IOobject io
+    (
+        lduAssemblyName_,
+        psi_.mesh().time().timeName(),
+        psi_.mesh(),
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    );
+
+    UPtrList<lduMesh> uMeshPtr(nMatrices());
+
+    UPtrList<GeometricField<Type, fvPatchField, volMesh>>
+        uFieldPtr(nMatrices());
+
+    for (label fieldi = 0; fieldi < nMatrices(); fieldi++)
+    {
+        const fvMesh& meshi = this->psi(fieldi).mesh();
+        uMeshPtr.set
+        (
+            fieldi,
+            &const_cast<fvMesh&>(meshi)
+        );
+        uFieldPtr.set(fieldi, &this->psi(fieldi));
+    }
+
+    if (!ptr)
+    {
+        lduPrimitiveMeshAssembly* lduAssemMeshPtr =
+            new lduPrimitiveMeshAssembly(io, uMeshPtr);
+
+        lduAssemMeshPtr->store();
+        lduAssemMeshPtr->update(uFieldPtr);
+
+        Info
+            << "Creating lduPrimitiveAssembly: " << lduAssemblyName_ << endl;
+    }
+    else if
+    (
+        psi_.mesh().changing() && !psi_.mesh().upToDatePoints(*ptr)
+    )
+    {
+        // Clear losortPtr_, ownerStartPtr_, losortStartPtr_
+        ptr->lduAddr().clearOut();
+        ptr->update(uFieldPtr);
+        psi_.mesh().setUpToDatePoints(*ptr);
+
+        Info
+            << "Updating lduPrimitiveAssembly: " << lduAssemblyName_ << endl;
+    }
+    else
+    {
+        Info
+            << "Using lduPrimitiveAssembly: " << lduAssemblyName_ << endl;
+    }
+}
+
+
 template<class Type>
 void Foam::fvMatrix<Type>::setValues
 (
@@ -555,6 +1141,36 @@ void Foam::fvMatrix<Type>::setReferences
 }
 
 
+template<class Type>
+void Foam::fvMatrix<Type>::addFvMatrix(fvMatrix& matrix)
+{
+    subMatrices_.append(matrix.clone());
+    ++nMatrix_;
+
+    if (dimensions_ != matrix.dimensions())
+    {
+        FatalErrorInFunction
+            << "incompatible dimensions for matrix addition "
+            << endl << "    "
+            << "[" << dimensions_  << " ] "
+            << " [" << matrix.dimensions() << " ]"
+            << abort(FatalError);
+    }
+
+    for (label fieldi = 0; fieldi < nMatrices(); fieldi++)
+    {
+        label id = checkImplicit(fieldi);
+        if (id > 0)
+        {
+            break;
+        }
+    }
+
+    internalCoeffs_.setSize(0);
+    boundaryCoeffs_.setSize(0);
+}
+
+
 template<class Type>
 void Foam::fvMatrix<Type>::relax(const scalar alpha)
 {
@@ -903,6 +1519,7 @@ Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Type>::H1() const
 }
 
 
+
 template<class Type>
 Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh>>
 Foam::fvMatrix<Type>::
@@ -917,6 +1534,14 @@ flux() const
             << abort(FatalError);
     }
 
+    if (nMatrices() > 1)
+    {
+        FatalErrorInFunction
+            << "Flux requested but " << psi_.name()
+            << " can't handle multiple fvMatrix."
+            << abort(FatalError);
+    }
+
     // construct GeometricField<Type, fvsPatchField, surfaceMesh>
     tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> tfieldFlux
     (
@@ -950,30 +1575,49 @@ flux() const
 
     FieldField<Field, Type> InternalContrib = internalCoeffs_;
 
-    forAll(InternalContrib, patchi)
+    label fieldi = 0;
+    if (!useImplicit_)
     {
-        InternalContrib[patchi] =
-            cmptMultiply
-            (
-                InternalContrib[patchi],
-                psi_.boundaryField()[patchi].patchInternalField()
-            );
+        forAll(InternalContrib, patchi)
+        {
+            InternalContrib[patchi] =
+                cmptMultiply
+                (
+                    InternalContrib[patchi],
+                    psi_.boundaryField()[patchi].patchInternalField()
+                );
+        }
+    }
+    else
+    {
+        FieldField<Field, Type> fluxInternalContrib(internalCoeffs_);
+
+        mapContributions(fieldi, fluxInternalContrib, InternalContrib, true);
     }
 
     FieldField<Field, Type> NeighbourContrib = boundaryCoeffs_;
 
-    forAll(NeighbourContrib, patchi)
+    if (!useImplicit_)
     {
-        if (psi_.boundaryField()[patchi].coupled())
+        forAll(NeighbourContrib, patchi)
         {
-            NeighbourContrib[patchi] =
-                cmptMultiply
-                (
-                    NeighbourContrib[patchi],
-                    psi_.boundaryField()[patchi].patchNeighbourField()
-                );
+            if (psi_.boundaryField()[patchi].coupled())
+            {
+                NeighbourContrib[patchi] =
+                    cmptMultiply
+                    (
+                        NeighbourContrib[patchi],
+                        psi_.boundaryField()[patchi].patchNeighbourField()
+                    );
+            }
         }
     }
+    else
+    {
+        FieldField<Field, Type> fluxBoundaryContrib(boundaryCoeffs_);
+
+        mapContributions(fieldi, fluxBoundaryContrib, NeighbourContrib, false);
+    }
 
     typename GeometricField<Type, fvsPatchField, surfaceMesh>::
         Boundary& ffbf = fieldFlux.boundaryFieldRef();
@@ -981,6 +1625,7 @@ flux() const
     forAll(ffbf, patchi)
     {
         ffbf[patchi] = InternalContrib[patchi] - NeighbourContrib[patchi];
+        //DebugVar(gSum(ffbf[patchi]))
     }
 
     if (faceFluxCorrectionPtr_)
@@ -1039,6 +1684,9 @@ void Foam::fvMatrix<Type>::operator=(const fvMatrix<Type>& fvmv)
             new GeometricField<Type, fvsPatchField, surfaceMesh>
         (*fvmv.faceFluxCorrectionPtr_);
     }
+
+    useImplicit_ = fvmv.useImplicit_;
+    lduAssemblyName_ = fvmv.lduAssemblyName_;
 }
 
 
@@ -1076,6 +1724,10 @@ void Foam::fvMatrix<Type>::operator+=(const fvMatrix<Type>& fvmv)
     internalCoeffs_ += fvmv.internalCoeffs_;
     boundaryCoeffs_ += fvmv.boundaryCoeffs_;
 
+    useImplicit_ = fvmv.useImplicit_;
+    lduAssemblyName_ = fvmv.lduAssemblyName_;
+    nMatrix_ = fvmv.nMatrix_;
+
     if (faceFluxCorrectionPtr_ && fvmv.faceFluxCorrectionPtr_)
     {
         *faceFluxCorrectionPtr_ += *fvmv.faceFluxCorrectionPtr_;
@@ -1110,6 +1762,10 @@ void Foam::fvMatrix<Type>::operator-=(const fvMatrix<Type>& fvmv)
     internalCoeffs_ -= fvmv.internalCoeffs_;
     boundaryCoeffs_ -= fvmv.boundaryCoeffs_;
 
+    useImplicit_ = fvmv.useImplicit_;
+    lduAssemblyName_ = fvmv.lduAssemblyName_;
+    nMatrix_ = fvmv.nMatrix_;
+
     if (faceFluxCorrectionPtr_ && fvmv.faceFluxCorrectionPtr_)
     {
         *faceFluxCorrectionPtr_ -= *fvmv.faceFluxCorrectionPtr_;

src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -51,6 +51,8 @@ SourceFiles
 #include "dimensionedTypes.H"
 #include "zero.H"
 #include "className.H"
+#include "lduPrimitiveMeshAssembly.H"
+#include "lduMesh.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -108,6 +110,7 @@ tmp<GeometricField<Type, fvPatchField, volMesh>> operator&
 );
 
 
+
 /*---------------------------------------------------------------------------*\
                            Class fvMatrix Declaration
 \*---------------------------------------------------------------------------*/
@@ -124,6 +127,18 @@ class fvMatrix
         //  Converted into a non-const reference at the point of solution.
         const GeometricField<Type, fvPatchField, volMesh>& psi_;
 
+        //- Originating fvMatrices when assembling matrices. Empty if not used.
+        PtrList<fvMatrix<Type>> subMatrices_;
+
+        //- Is the fvMatrix using implicit formulation
+        bool useImplicit_;
+
+        //- Name of the lduAssembly
+        word lduAssemblyName_;
+
+        //- Number of fvMatrices added to this
+        label nMatrix_;
+
         //- Dimension set
         dimensionSet dimensions_;
 
@@ -148,6 +163,7 @@ protected:
     //- Declare friendship with the fvSolver class
     friend class fvSolver;
 
+
     // Protected Member Functions
 
         //- Add patch contribution to internal field
@@ -185,6 +201,12 @@ protected:
         ) const;
 
 
+        // Implicit helper functions
+
+            //- Name the implicit assembly addressing
+            label checkImplicit(const label fieldI = 0);
+
+
         // Matrix completion functionality
 
             void addBoundaryDiag
@@ -285,10 +307,121 @@ public:
 
         // Access
 
-            const GeometricField<Type, fvPatchField, volMesh>& psi() const
-            {
-                return psi_;
-            }
+            // Coupling
+
+                label nMatrices() const
+                {
+                    return (nMatrix_ == 0 ? 1 : nMatrix_);
+                }
+
+                const fvMatrix<Type>& matrix(const label i) const
+                {
+                    return (nMatrix_ == 0 ? *this : subMatrices_[i]);
+                }
+
+                fvMatrix<Type>& matrix(const label i)
+                {
+                    return (nMatrix_ == 0 ? *this : subMatrices_[i]);
+                }
+
+                label globalPatchID
+                (
+                    const label fieldi,
+                    const label patchi
+                ) const
+                {
+                    if (!lduMeshPtr())
+                    {
+                        return patchi;
+                    }
+                    else
+                    {
+                        return lduMeshPtr()->patchMap()[fieldi][patchi];
+                    }
+                }
+
+                //- Transfer lower, upper, diag and source to this fvMatrix
+                void transferFvMatrixCoeffs();
+
+                //- Create or update ldu assembly
+                void createOrUpdateLduPrimitiveAssembly();
+
+                //- Access to lduPrimitiveMeshAssembly
+                lduPrimitiveMeshAssembly* lduMeshPtr();
+
+                //- Const Access to lduPrimitiveMeshAssembly
+                const lduPrimitiveMeshAssembly* lduMeshPtr() const;
+
+                //- Manipulate matrix
+                void manipulateMatrix(direction cmp);
+
+                //- Manipulate boundary/internal coeffs for coupling
+                void setBounAndInterCoeffs();
+
+                //- Set interfaces
+                void setInterfaces
+                (
+                    lduInterfaceFieldPtrsList&,
+                    PtrDynList<lduInterfaceField>& newInterfaces
+                );
+
+                //- Add internal and boundary contribution to local patches
+                void mapContributions
+                (
+                    label fieldi,
+                    const FieldField<Field, Type>& fluxContrib,
+                    FieldField<Field, Type>& contrib,
+                    bool internal
+                ) const;
+
+                //- Return optional lduAdressing
+                const lduPrimitiveMeshAssembly& lduMeshAssembly()
+                {
+                    return *lduMeshPtr();
+                }
+
+                //- Return psi
+                const GeometricField<Type, fvPatchField, volMesh>& psi
+                (
+                    const label i = 0
+                ) const
+                {
+                    return
+                    (
+                        (i == 0 && nMatrix_ == 0) ? psi_ : matrix(i).psi()
+                    );
+                }
+
+
+                GeometricField<Type, fvPatchField, volMesh>& psi
+                (
+                    const label i = 0
+                )
+                {
+                    return
+                    (
+                        (i == 0 && nMatrix_ == 0)
+                      ? const_cast
+                        <
+                            GeometricField<Type, fvPatchField, volMesh>&
+                        >(psi_)
+                      : const_cast
+                        <
+                            GeometricField<Type, fvPatchField, volMesh>&
+                        >
+                        (
+                            matrix(i).psi()
+                        )
+                    );
+                }
+
+                //- Clear multiple fvMatrices
+                void clear()
+                {
+                    subMatrices_.clear();
+                    nMatrix_ = 0;
+                }
+
 
             const dimensionSet& dimensions() const
             {
@@ -333,7 +466,6 @@ public:
                 return boundaryCoeffs_;
             }
 
-
             //- Declare return type of the faceFluxCorrectionPtr() function
             typedef GeometricField<Type, fvsPatchField, surfaceMesh>
                 *surfaceTypeFieldPtr;
@@ -405,6 +537,9 @@ public:
                 const scalar value
             );
 
+            //- Add fvMatrix
+            void addFvMatrix(fvMatrix<Type>& matrix);
+
             //- Relax matrix (for steady-state solution).
             //  alpha = 1 : diagonally equal
             //  alpha < 1 : diagonally dominant

src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C

@@ -116,6 +116,12 @@ Foam::SolverPerformance<Type> Foam::fvMatrix<Type>::solveSegregated
     const dictionary& solverControls
 )
 {
+    if (useImplicit_)
+    {
+        FatalErrorInFunction
+            << "Implicit option is not allowed for type: " << Type::typeName
+            << exit(FatalError);
+    }
     if (debug)
     {
         Info.masterStream(this->mesh().comm())

src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C

@@ -30,6 +30,9 @@ License
 #include "extrapolatedCalculatedFvPatchFields.H"
 #include "profiling.H"
 #include "PrecisionAdaptor.H"
+#include "jumpCyclicFvPatchField.H"
+#include "cyclicPolyPatch.H"
+#include "cyclicAMIPolyPatch.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -82,6 +85,10 @@ Foam::fvMatrix<Foam::scalar>::solver
     scalarField saveDiag(diag());
     addBoundaryDiag(diag(), 0);
 
+    lduInterfaceFieldPtrsList interfaces =
+        psi_.boundaryField().scalarInterfaces();
+
+
     autoPtr<fvMatrix<scalar>::fvSolver> solverPtr
     (
         new fvMatrix<scalar>::fvSolver
@@ -93,7 +100,7 @@ Foam::fvMatrix<Foam::scalar>::solver
                 *this,
                 boundaryCoeffs_,
                 internalCoeffs_,
-                psi_.boundaryField().scalarInterfaces(),
+                interfaces,
                 solverControls
             )
         )
@@ -160,8 +167,29 @@ Foam::solverPerformance Foam::fvMatrix<Foam::scalar>::solveSegregated
             << endl;
     }
 
-    GeometricField<scalar, fvPatchField, volMesh>& psi =
-       const_cast<GeometricField<scalar, fvPatchField, volMesh>&>(psi_);
+    scalarField saveLower;
+    scalarField saveUpper;
+
+    if (useImplicit_)
+    {
+        createOrUpdateLduPrimitiveAssembly();
+
+        if (psi_.mesh().fluxRequired(psi_.name()))
+        {
+            // Save lower/upper for flux calculation
+            if (asymmetric())
+            {
+                saveLower = lower();
+            }
+            saveUpper = upper();
+        }
+
+        setLduMesh(*lduMeshPtr());
+        transferFvMatrixCoeffs();
+        setBounAndInterCoeffs();
+        direction cmpt = 0;
+        manipulateMatrix(cmpt);
+    }
 
     scalarField saveDiag(diag());
     addBoundaryDiag(diag(), 0);
@@ -169,16 +197,72 @@ Foam::solverPerformance Foam::fvMatrix<Foam::scalar>::solveSegregated
     scalarField totalSource(source_);
     addBoundarySource(totalSource, false);
 
+    lduInterfaceFieldPtrsList interfaces;
+    PtrDynList<lduInterfaceField> newInterfaces;
+    if (!useImplicit_)
+    {
+        interfaces = this->psi(0).boundaryField().scalarInterfaces();
+    }
+    else
+    {
+        setInterfaces(interfaces, newInterfaces);
+    }
+
+    tmp<scalarField> tpsi;
+    if (!useImplicit_)
+    {
+        GeometricField<scalar, fvPatchField, volMesh>& psi =
+           const_cast<GeometricField<scalar, fvPatchField, volMesh>&>(psi_);
+
+        tpsi = tmp<scalarField>(psi.primitiveFieldRef());
+    }
+    else
+    {
+        tpsi = tmp<scalarField>::New(lduAddr().size(), Zero);
+        scalarField& psi = tpsi.ref();
+
+        for (label fieldi = 0; fieldi < nMatrices(); fieldi++)
+        {
+            const label cellOffset = lduMeshPtr()->cellOffsets()[fieldi];
+            const auto& psiInternal = this->psi(fieldi).primitiveField();
+
+            forAll(psiInternal, localCellI)
+            {
+                psi[cellOffset + localCellI] = psiInternal[localCellI];
+            }
+        }
+    }
+    scalarField& psi = tpsi.constCast();
+
     // Solver call
     solverPerformance solverPerf = lduMatrix::solver::New
     (
-        psi.name(),
+        this->psi(0).name(),
         *this,
         boundaryCoeffs_,
         internalCoeffs_,
-        psi_.boundaryField().scalarInterfaces(),
+        interfaces,
         solverControls
-    )->solve(psi.primitiveFieldRef(), totalSource);
+    )->solve(psi, totalSource);
+
+    if (useImplicit_)
+    {
+        for (label fieldi = 0; fieldi < nMatrices(); fieldi++)
+        {
+            auto& psiInternal =
+                const_cast<GeometricField<scalar, fvPatchField, volMesh>&>
+                (
+                    this->psi(fieldi)
+                ).primitiveFieldRef();
+
+            const label cellOffset = lduMeshPtr()->cellOffsets()[fieldi];
+
+            forAll(psiInternal, localCellI)
+            {
+                psiInternal[localCellI] = psi[localCellI + cellOffset];
+            }
+        }
+    }
 
     if (solverPerformance::debug)
     {
@@ -187,9 +271,35 @@ Foam::solverPerformance Foam::fvMatrix<Foam::scalar>::solveSegregated
 
     diag() = saveDiag;
 
-    psi.correctBoundaryConditions();
+    if (useImplicit_)
+    {
+        if (psi_.mesh().fluxRequired(psi_.name()))
+        {
+            // Restore lower/upper
+            if (asymmetric())
+            {
+                lower().setSize(saveLower.size());
+                lower() = saveLower;
+            }
 
-    psi.mesh().setSolverPerformance(psi.name(), solverPerf);
+            upper().setSize(saveUpper.size());
+            upper() = saveUpper;
+        }
+        // Set the original lduMesh
+        setLduMesh(psi_.mesh());
+    }
+
+    for (label fieldi = 0; fieldi < nMatrices(); fieldi++)
+    {
+        auto& localPsi =
+            const_cast<GeometricField<scalar, fvPatchField, volMesh>&>
+            (
+                this->psi(fieldi)
+            );
+
+        localPsi.correctBoundaryConditions();
+        localPsi.mesh().setSolverPerformance(localPsi.name(), solverPerf);
+    }
 
     return solverPerf;
 }

src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.H

@@ -59,22 +59,13 @@ void fvMatrix<scalar>::setComponentReference
 );
 
 template<>
-autoPtr<fvMatrix<scalar>::fvSolver> fvMatrix<scalar>::solver
-(
-    const dictionary&
-);
+autoPtr<fvMatrix<scalar>::fvSolver> fvMatrix<scalar>::solver(const dictionary&);
 
 template<>
-solverPerformance fvMatrix<scalar>::fvSolver::solve
-(
-    const dictionary&
-);
+solverPerformance fvMatrix<scalar>::fvSolver::solve(const dictionary&);
 
 template<>
-solverPerformance fvMatrix<scalar>::solveSegregated
-(
-    const dictionary&
-);
+solverPerformance fvMatrix<scalar>::solveSegregated(const dictionary&);
 
 template<>
 tmp<scalarField> fvMatrix<scalar>::residual() const;

src/finiteVolume/fvMesh/fvMesh.C

@@ -697,12 +697,20 @@ const Foam::lduAddressing& Foam::fvMesh::lduAddr() const
             << nFaces() << endl;
 
         lduPtr_ = new fvMeshLduAddressing(*this);
+
+        return *lduPtr_;
     }
 
     return *lduPtr_;
 }
 
 
+Foam::lduInterfacePtrsList Foam::fvMesh::interfaces() const
+{
+    return boundary().interfaces();
+}
+
+
 void Foam::fvMesh::mapFields(const mapPolyMesh& meshMap)
 {
     DebugInFunction

src/finiteVolume/fvMesh/fvMesh.H

@@ -67,6 +67,8 @@ SourceFiles
 #include "className.H"
 #include "SolverPerformance.H"
 
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -309,10 +311,7 @@ public:
 
             //- Return a list of pointers for each patch
             //  with only those pointing to interfaces being set
-            virtual lduInterfacePtrsList interfaces() const
-            {
-                return boundary().interfaces();
-            }
+            virtual lduInterfacePtrsList interfaces() const;
 
             //- Return communicator used for parallel communication
             virtual label comm() const

src/finiteVolume/fvMesh/fvPatches/basic/coupled/coupledFvPatch.H

@@ -124,12 +124,21 @@ public:
         // Interface transfer functions
 
             //- Return the values of the given internal data adjacent to
-            //  the interface as a field
+            //- the interface as a field
             virtual tmp<labelField> interfaceInternalField
             (
                 const labelUList& internalData
             ) const = 0;
 
+
+            //- Return the values of the given internal data adjacent to
+            //- the interface as a field using a mapping faceCell
+            virtual tmp<labelField> interfaceInternalField
+            (
+                const labelUList& internalData,
+                const labelUList& faceCell
+            ) const = 0;
+
             //- Inherit initInternalFieldTransfer from lduInterface
             using lduInterface::initInternalFieldTransfer;
 
@@ -141,6 +150,15 @@ public:
             ) const
             {}
 
+            //- Initialise neighbour field transfer using faceCells map
+            virtual void initInternalFieldTransfer
+            (
+                const Pstream::commsTypes commsType,
+                labelUList& iF,
+                const labelUList& faceCell
+            ) const
+            {}
+
             //- Return neighbour field
             virtual tmp<labelField> internalFieldTransfer
             (

src/finiteVolume/fvMesh/fvPatches/constraint/cyclic/cyclicFvPatch.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -101,6 +102,16 @@ Foam::tmp<Foam::labelField> Foam::cyclicFvPatch::interfaceInternalField
 }
 
 
+Foam::tmp<Foam::labelField> Foam::cyclicFvPatch::interfaceInternalField
+(
+    const labelUList& internalData,
+    const labelUList& faceCells
+) const
+{
+    return patchInternalField(internalData, faceCells);
+}
+
+
 Foam::tmp<Foam::labelField> Foam::cyclicFvPatch::internalFieldTransfer
 (
     const Pstream::commsTypes commsType,

src/finiteVolume/fvMesh/fvPatches/constraint/cyclic/cyclicFvPatch.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -87,6 +88,7 @@ public:
 
     // Member functions
 
+
         // Access
 
             //- Return local reference cast into the cyclic patch
@@ -149,18 +151,29 @@ public:
         // Interface transfer functions
 
             //- Return the values of the given internal data adjacent to
-            //  the interface as a field
+            //- the interface as a field
             virtual tmp<labelField> interfaceInternalField
             (
                 const labelUList& internalData
             ) const;
 
+
+            //- Return the values of the given internal data adjacent to
+            //- the interface as a field using a mapping faceCell
+            virtual tmp<labelField> interfaceInternalField
+            (
+                const labelUList& internalData,
+                const labelUList& faceCell
+            ) const;
+
+
             //- Return neighbour field
             virtual tmp<labelField> internalFieldTransfer
             (
                 const Pstream::commsTypes commsType,
                 const labelUList& internalData
             ) const;
+
 };
 
 

src/finiteVolume/fvMesh/fvPatches/constraint/cyclicACMI/cyclicACMIFvPatch.C

@@ -168,6 +168,68 @@ Foam::cyclicACMIFvPatch::cyclicACMIFvPatch
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+// void Foam::cyclicACMIFvPatch::newInternalProcFaces
+// (
+//     label& newFaces,
+//     label& newProcFaces
+// ) const
+// {
+//     const List<labelList>& addSourceFaces =
+//         cyclicACMIPolyPatch_.AMI().srcAddress();
+//
+//     const scalarField& fMask = cyclicACMIPolyPatch_.srcMask();
+//
+//     // Add new faces as many weights for AMI
+//     forAll (addSourceFaces, faceI)
+//     {
+//         if (fMask[faceI] > cyclicACMIPolyPatch_.tolerance_)
+//         {
+//             const labelList& nbrFaceIs = addSourceFaces[faceI];
+//
+//             forAll (nbrFaceIs, j)
+//             {
+//                 label nbrFaceI = nbrFaceIs[j];
+//
+//                 if (nbrFaceI < neighbPatch().size())
+//                 {
+//                     // local faces
+//                     newFaces++;
+//                 }
+//                 else
+//                 {
+//                     // Proc faces
+//                     newProcFaces++;
+//                 }
+//             }
+//         }
+//     }
+// }
+
+
+// Foam::refPtr<Foam::labelListList>
+// Foam::cyclicACMIFvPatch::mapCollocatedFaces() const
+// {
+//     const scalarField& fMask = cyclicACMIPolyPatch_.srcMask();
+//     const labelListList& srcFaces = cyclicACMIPolyPatch_.AMI().srcAddress();
+//     labelListList dOverFaces;
+//
+//     dOverFaces.setSize(srcFaces.size());
+//     forAll (dOverFaces, faceI)
+//     {
+//         if (fMask[faceI] > cyclicACMIPolyPatch_.tolerance_)
+//         {
+//             dOverFaces[faceI].setSize(srcFaces[faceI].size());
+//
+//             forAll (dOverFaces[faceI], subFaceI)
+//             {
+//                 dOverFaces[faceI][subFaceI] = srcFaces[faceI][subFaceI];
+//             }
+//         }
+//     }
+//     return refPtr<labelListList>(new labelListList(dOverFaces));
+// }
+
+
 bool Foam::cyclicACMIFvPatch::coupled() const
 {
     return Pstream::parRun() || (this->size() && neighbFvPatch().size());
@@ -227,6 +289,16 @@ Foam::tmp<Foam::labelField> Foam::cyclicACMIFvPatch::interfaceInternalField
 }
 
 
+Foam::tmp<Foam::labelField> Foam::cyclicACMIFvPatch::interfaceInternalField
+(
+    const labelUList& internalData,
+    const labelUList& faceCells
+) const
+{
+    return patchInternalField(internalData, faceCells);
+}
+
+
 Foam::tmp<Foam::labelField> Foam::cyclicACMIFvPatch::internalFieldTransfer
 (
     const Pstream::commsTypes commsType,
@@ -344,5 +416,4 @@ void Foam::cyclicACMIFvPatch::movePoints()
     }
 }
 
-
 // ************************************************************************* //

src/finiteVolume/fvMesh/fvPatches/constraint/cyclicACMI/cyclicACMIFvPatch.H

@@ -97,6 +97,15 @@ public:
 
     // Member functions
 
+//         // Implicit treatment functions
+//
+//             //- Return number of new internal sub-faces and new proc faces
+//             virtual void newInternalProcFaces(label&, label&) const;
+//
+//             //- Return collocated faces
+//             virtual refPtr<labelListList> mapCollocatedFaces() const;
+
+
         // Access
 
             //- Return local reference cast into the cyclic patch
@@ -169,9 +178,9 @@ public:
                 );
             }
 
-            //- Return true if this patch is coupled. This is equivalent
-            //- to the coupledPolyPatch::coupled() if parallel running or
-            //- both sides present, false otherwise
+            //- Return true if this patch is coupled.
+            //  This is equivalent to coupledPolyPatch::coupled()
+            //  if parallel running or both sides present, false otherwise
             virtual bool coupled() const;
 
             //- Return delta (P to N) vectors across coupled patch
@@ -214,12 +223,21 @@ public:
                 const labelUList& internalData
             ) const;
 
+            //- Return the values of the given internal data adjacent to
+            //- the interface as a field using a mapping faceCell
+            virtual tmp<labelField> interfaceInternalField
+            (
+                const labelUList& internalData,
+                const labelUList& faceCells
+            ) const;
+
             //- Return neighbour field
             virtual tmp<labelField> internalFieldTransfer
             (
                 const Pstream::commsTypes commsType,
                 const labelUList& internalData
             ) const;
+
 };
 
 

src/finiteVolume/fvMesh/fvPatches/constraint/cyclicAMI/cyclicAMIFvPatch.C

@@ -51,6 +51,38 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+// void Foam::cyclicAMIFvPatch::newInternalProcFaces
+// (
+//     label& newFaces,
+//     label& newProcFaces
+// ) const
+// {
+//     const labelListList& addSourceFaces = AMI().srcAddress();
+//
+//     // Add new faces as many weights for AMI
+//     forAll (addSourceFaces, faceI)
+//     {
+//         const labelList& nbrFaceIs = addSourceFaces[faceI];
+//
+//         forAll (nbrFaceIs, j)
+//         {
+//             label nbrFaceI = nbrFaceIs[j];
+//
+//             if (nbrFaceI < neighbPatch().size())
+//             {
+//                 // local faces
+//                 newFaces++;
+//             }
+//             else
+//             {
+//                 // Proc faces
+//                 newProcFaces++;
+//             }
+//         }
+//     }
+// }
+
+
 bool Foam::cyclicAMIFvPatch::coupled() const
 {
     return
@@ -190,6 +222,16 @@ Foam::tmp<Foam::labelField> Foam::cyclicAMIFvPatch::interfaceInternalField
 }
 
 
+Foam::tmp<Foam::labelField> Foam::cyclicAMIFvPatch::interfaceInternalField
+(
+    const labelUList& internalData,
+    const labelUList& faceCells
+) const
+{
+    return patchInternalField(internalData, faceCells);
+}
+
+
 Foam::tmp<Foam::labelField> Foam::cyclicAMIFvPatch::internalFieldTransfer
 (
     const Pstream::commsTypes commsType,
@@ -270,5 +312,4 @@ void Foam::cyclicAMIFvPatch::movePoints()
     }
 }
 
-
 // ************************************************************************* //

src/finiteVolume/fvMesh/fvPatches/constraint/cyclicAMI/cyclicAMIFvPatch.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -102,6 +103,41 @@ public:
 
         // Access
 
+//             // Implicit treatment functions
+//
+//                 //- Return number of new internal of this polyPatch faces
+//                 virtual void newInternalProcFaces(label&, label&) const;
+//
+//                 //- Return nbrCells
+//                 virtual const labelUList& nbrCells() const
+//                 {
+//                     return cyclicAMIPolyPatch_.neighbPatch().faceCells();
+//                 }
+//
+//                 //- Return nbr patch ID
+//                 virtual label neighbPolyPatchID() const
+//                 {
+//                     return neighbPatchID();
+//                 }
+//
+//                 //- Return collocated faces map
+//                 virtual refPtr<labelListList> mapCollocatedFaces() const
+//                 {
+//                     const labelListList& sourceFaces =
+//                         cyclicAMIPolyPatch_.AMI().srcAddress();
+//                     return refPtr<labelListList>
+//                     (
+//                         new labelListList(sourceFaces)
+//                     );
+//                 }
+//
+//                 //- Return implicit master
+//                 virtual bool masterImplicit() const
+//                 {
+//                     return owner();
+//                 }
+
+
             //- Return local reference cast into the cyclic patch
             const cyclicAMIPolyPatch& cyclicAMIPatch() const
             {
@@ -167,6 +203,7 @@ public:
                 );
             }
 
+
             //- Return true if this patch is coupled. This is equivalent
             //- to the coupledPolyPatch::coupled() if parallel running or
             //- both sides present, false otherwise
@@ -205,12 +242,23 @@ public:
                 const labelUList& internalData
             ) const;
 
+
+            //- Return the values of the given internal data adjacent to
+            //- the interface as a field using a mapping faceCell
+            virtual tmp<labelField> interfaceInternalField
+            (
+                const labelUList& internalData,
+                const labelUList& faceCells
+            ) const;
+
+
             //- Return neighbour field
             virtual tmp<labelField> internalFieldTransfer
             (
                 const Pstream::commsTypes commsType,
                 const labelUList& internalData
             ) const;
+
 };
 
 

src/finiteVolume/fvMesh/fvPatches/constraint/processor/processorFvPatch.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2014 OpenFOAM Foundation
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -110,6 +111,16 @@ Foam::tmp<Foam::labelField> Foam::processorFvPatch::interfaceInternalField
 }
 
 
+Foam::tmp<Foam::labelField> Foam::processorFvPatch::interfaceInternalField
+(
+    const labelUList& internalData,
+    const labelUList& faceCells
+) const
+{
+    return patchInternalField(internalData, faceCells);
+}
+
+
 void Foam::processorFvPatch::initInternalFieldTransfer
 (
     const Pstream::commsTypes commsType,
@@ -120,6 +131,17 @@ void Foam::processorFvPatch::initInternalFieldTransfer
 }
 
 
+void Foam::processorFvPatch::initInternalFieldTransfer
+(
+    const Pstream::commsTypes commsType,
+    const labelUList& iF,
+    const labelUList& faceCells
+) const
+{
+    send(commsType, interfaceInternalField(iF, faceCells)());
+}
+
+
 Foam::tmp<Foam::labelField> Foam::processorFvPatch::internalFieldTransfer
 (
     const Pstream::commsTypes commsType,