From c8538ee49ec16aa20277406e26a4bcea149d021e Mon Sep 17 00:00:00 2001
From: sergio <s.ferraris@opencfd.co.uk>
Date: Mon, 31 Jan 2022 12:10:08 -0800
Subject: [PATCH] ENH: Adding iterative solver to S2S equations

ENH: Adding report after matrix smoothing
---
 .../preProcessing/viewFactorsGen/shootRays.H  |   4 +-
 .../viewFactorsGen/viewFactorsGen.C           |  17 +-
 .../lduCalculatedProcessorField.C             | 236 ++++++
 .../lduCalculatedProcessorField.H             | 254 ++++++
 .../lduPrimitiveProcessorInterface.C          |  16 +-
 .../lduPrimitiveProcessorInterface.H          |  12 +-
 .../radiationModels/viewFactor/viewFactor.C   | 736 ++++++++++++++----
 .../radiationModels/viewFactor/viewFactor.H   |  39 +
 8 files changed, 1146 insertions(+), 168 deletions(-)
 create mode 100644 src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduCalculatedProcessorField/lduCalculatedProcessorField.C
 create mode 100644 src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduCalculatedProcessorField/lduCalculatedProcessorField.H

diff --git a/applications/utilities/preProcessing/viewFactorsGen/shootRays.H b/applications/utilities/preProcessing/viewFactorsGen/shootRays.H
index 9f751ecd20..1df6e108f0 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/shootRays.H
+++ b/applications/utilities/preProcessing/viewFactorsGen/shootRays.H
@@ -4,7 +4,7 @@
 // Maximum length for dynamicList
 const label maxDynListLength
 (
-    viewFactorDict.getOrDefault<label>("maxDynListLength", 100000)
+    viewFactorDict.getOrDefault<label>("maxDynListLength", 1000000)
 );
 
 for (const int proci : Pstream::allProcs())
@@ -47,7 +47,7 @@ for (const int proci : Pstream::allProcs())
                     const vector& remA = remoteArea[j];
                     const label& remAgg = remoteAgg[j];
 
-                    const vector& d = remFc - fc;
+                    const vector d(remFc - fc);
 
                     if (((d & fA) < 0.) && ((d & remA) > 0))
                     {
diff --git a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
index a32b878ec9..fafdfb811b 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
+++ b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
@@ -408,7 +408,7 @@ int main(int argc, char *argv[])
                 List<point> availablePoints
                 (
                     upp.faceCentres().size()
-                + upp.localPoints().size()
+                  + upp.localPoints().size()
                 );
 
                 SubList<point>
@@ -863,21 +863,6 @@ int main(int argc, char *argv[])
     IOglobalFaceFaces.write();
 
 
-    labelListIOList IOvisibleFaceFaces
-    (
-        IOobject
-        (
-            "visibleFaceFaces",
-            mesh.facesInstance(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        std::move(visibleFaceFaces)
-    );
-    IOvisibleFaceFaces.write();
-
     IOmapDistribute IOmapDist
     (
         IOobject
diff --git a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduCalculatedProcessorField/lduCalculatedProcessorField.C b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduCalculatedProcessorField/lduCalculatedProcessorField.C
new file mode 100644
index 0000000000..1108e24a4c
--- /dev/null
+++ b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduCalculatedProcessorField/lduCalculatedProcessorField.C
@@ -0,0 +1,236 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "lduCalculatedProcessorField.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::lduCalculatedProcessorField<Type>::lduCalculatedProcessorField
+(
+    const lduInterface& interface,
+    const Field<Type>& iF
+)
+:
+    LduInterfaceField<Type>(interface),
+    procInterface_(refCast<const lduPrimitiveProcessorInterface>(interface)),
+    field_(iF),
+    sendBuf_(procInterface_.faceCells().size()),
+    receiveBuf_(procInterface_.faceCells().size()),
+    scalarSendBuf_(procInterface_.faceCells().size()),
+    scalarReceiveBuf_(procInterface_.faceCells().size()),
+    outstandingSendRequest_(-1),
+    outstandingRecvRequest_(-1)
+{}
+
+
+template<class Type>
+Foam::lduCalculatedProcessorField<Type>::lduCalculatedProcessorField
+(
+    const lduCalculatedProcessorField<Type>& ptf
+)
+:
+    LduInterfaceField<Type>(refCast<const lduInterface>(ptf)),
+    procInterface_(ptf.procInterface_),
+    field_(ptf.field_),
+    sendBuf_(procInterface_.faceCells().size()),
+    receiveBuf_(procInterface_.faceCells().size()),
+    scalarSendBuf_(procInterface_.faceCells().size()),
+    scalarReceiveBuf_(procInterface_.faceCells().size()),
+    outstandingSendRequest_(-1),
+    outstandingRecvRequest_(-1)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+bool Foam::lduCalculatedProcessorField<Type>::ready() const
+{
+    if
+    (
+        this->outstandingSendRequest_ >= 0
+     && this->outstandingSendRequest_ < Pstream::nRequests()
+    )
+    {
+        if (!UPstream::finishedRequest(this->outstandingSendRequest_))
+        {
+            return false;
+        }
+    }
+    this->outstandingSendRequest_ = -1;
+
+    if
+    (
+        this->outstandingRecvRequest_ >= 0
+     && this->outstandingRecvRequest_ < Pstream::nRequests()
+    )
+    {
+        if (!UPstream::finishedRequest(this->outstandingRecvRequest_))
+        {
+            return false;
+        }
+    }
+    this->outstandingRecvRequest_ = -1;
+
+    return true;
+}
+
+
+template<class Type>
+void Foam::lduCalculatedProcessorField<Type>::initInterfaceMatrixUpdate
+(
+    solveScalarField& result,
+    const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
+    const solveScalarField& psiInternal,
+    const scalarField& coeffs,
+    const direction cmpt,
+    const Pstream::commsTypes commsType
+) const
+{
+    // Bypass patchInternalField since uses fvPatch addressing
+    const labelList& fc = lduAddr.patchAddr(patchId);
+
+    scalarSendBuf_.setSize(fc.size());
+    forAll(fc, i)
+    {
+        scalarSendBuf_[i] = psiInternal[fc[i]];
+    }
+
+    if (!this->ready())
+    {
+        FatalErrorInFunction
+            << "On patch "
+            << " outstanding request."
+            << abort(FatalError);
+    }
+
+
+    scalarReceiveBuf_.setSize(scalarSendBuf_.size());
+    outstandingRecvRequest_ = UPstream::nRequests();
+
+    UIPstream::read
+    (
+        Pstream::commsTypes::nonBlocking,
+        procInterface_.neighbProcNo(),
+        scalarReceiveBuf_.data_bytes(),
+        scalarReceiveBuf_.size_bytes(),
+        procInterface_.tag(),
+        procInterface_.comm()
+    );
+
+    outstandingSendRequest_ = UPstream::nRequests();
+
+    UOPstream::write
+    (
+        Pstream::commsTypes::nonBlocking,
+        procInterface_.neighbProcNo(),
+        scalarSendBuf_.cdata_bytes(),
+        scalarSendBuf_.size_bytes(),
+        procInterface_.tag(),
+        procInterface_.comm()
+    );
+
+    const_cast<lduInterfaceField&>
+    (
+        static_cast<const lduInterfaceField&>(*this)
+    ).updatedMatrix() = false;
+}
+
+
+template<class Type>
+void Foam::lduCalculatedProcessorField<Type>::addToInternalField
+(
+    solveScalarField& result,
+    const bool add,
+    const scalarField& coeffs,
+    const solveScalarField& vals
+) const
+{
+    const labelUList& faceCells = this->procInterface_.faceCells();
+
+    if (add)
+    {
+        forAll(faceCells, elemI)
+        {
+            result[faceCells[elemI]] += coeffs[elemI]*vals[elemI];
+        }
+    }
+    else
+    {
+        forAll(faceCells, elemI)
+        {
+            result[faceCells[elemI]] -= coeffs[elemI]*vals[elemI];
+        }
+    }
+}
+
+
+template<class Type>
+void Foam::lduCalculatedProcessorField<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 commsType
+) const
+{
+    if (this->updatedMatrix())
+    {
+        return;
+    }
+
+    if
+    (
+        outstandingRecvRequest_ >= 0
+     && outstandingRecvRequest_ < Pstream::nRequests()
+    )
+    {
+        UPstream::waitRequest(outstandingRecvRequest_);
+    }
+    // Recv finished so assume sending finished as well.
+    outstandingSendRequest_ = -1;
+    outstandingRecvRequest_ = -1;
+
+    // Consume straight from scalarReceiveBuf_. Note use of our own
+    // helper to avoid using fvPatch addressing
+    addToInternalField(result, !add, coeffs, scalarReceiveBuf_);
+
+    const_cast<lduInterfaceField&>
+    (
+        static_cast<const lduInterfaceField&>(*this)
+    ).updatedMatrix() = true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduCalculatedProcessorField/lduCalculatedProcessorField.H b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduCalculatedProcessorField/lduCalculatedProcessorField.H
new file mode 100644
index 0000000000..f9a36b1bfb
--- /dev/null
+++ b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduCalculatedProcessorField/lduCalculatedProcessorField.H
@@ -0,0 +1,254 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::lduCalculatedProcessorField
+
+Group
+    grpGenericBoundaryConditions
+
+Description
+    A lduProcessorField type bypassing coupledFvPatchField and holding
+    a reference to the Field<Type>.
+
+    Used to add updateInterfaceMatrix capabilities to a lduMatrix
+    which is fully uncoupled from the fvMesh.
+
+    Its functionality is purely to init and update the processor interfaces.
+
+SourceFiles
+    lduCalculatedProcessorField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef lduCalculatedProcessorField_H
+#define lduCalculatedProcessorField_H
+
+#include "lduPrimitiveProcessorInterface.H"
+#include "processorLduInterfaceField.H"
+#include "LduInterfaceField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+               Class lduCalculatedProcessorField Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class lduCalculatedProcessorField
+:
+    public LduInterfaceField<Type>,
+    public processorLduInterfaceField
+{
+protected:
+
+    // Protected Data
+
+        //- Local reference cast into the interface
+         const lduPrimitiveProcessorInterface& procInterface_;
+
+        //- Local Field
+         const Field<Type>& field_;
+
+        // Sending and receiving
+
+            //- Send buffer
+            mutable Field<Type> sendBuf_;
+
+            //- Receive buffer
+            mutable Field<Type> receiveBuf_;
+
+            //- Scalar send buffer
+            mutable solveScalarField scalarSendBuf_;
+
+            //- Scalar receive buffer
+            mutable solveScalarField scalarReceiveBuf_;
+
+            //- Outstanding request
+            mutable label outstandingSendRequest_;
+
+            //- Outstanding request
+            mutable label outstandingRecvRequest_;
+
+
+    // Protected Member Functions
+
+        void addToInternalField
+        (
+            solveScalarField& result,
+            const bool add,
+            const scalarField& coeffs,
+            const solveScalarField& vals
+        ) const;
+
+
+public:
+
+    //- Runtime type information
+    ClassName("lduCalculatedProcessorField");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        lduCalculatedProcessorField
+        (
+            const lduInterface& interface,
+            const Field<Type>&
+        );
+
+        //- Construct as copy
+        lduCalculatedProcessorField
+        (
+            const lduCalculatedProcessorField<Type>&
+        );
+
+
+    //- Destructor
+    virtual ~lduCalculatedProcessorField() = default;
+
+
+    // Member Functions
+
+    // Access
+
+        //- Return communicator used for comms
+        virtual label comm() const
+        {
+            return procInterface_.comm();
+        }
+
+        //- Return processor number
+        virtual int myProcNo() const
+        {
+            return procInterface_.myProcNo();
+        }
+
+        //- Return neighbour processor number
+        virtual int neighbProcNo() const
+        {
+            return procInterface_.myProcNo();
+        }
+
+        //- Is the transform required
+        virtual bool doTransform() const
+        {
+            return false;
+        }
+
+        //- Return face transformation tensor
+        virtual const tensorField& forwardT() const
+        {
+            return procInterface_.forwardT();
+        }
+
+        //- Return rank of component for transform
+        virtual int rank() const
+        {
+            return pTraits<Type>::rank;
+        }
+
+
+    // Evaluation
+
+        //- Is all data available
+        virtual bool ready() const;
+
+        //- Initialise neighbour matrix update
+        virtual void initInterfaceMatrixUpdate
+        (
+            solveScalarField& result,
+            const bool add,
+            const lduAddressing& lduAddr,
+            const label patchId,
+            const solveScalarField& psiInternal,
+            const scalarField& coeffs,
+            const direction cmpt,
+            const Pstream::commsTypes commsType
+        ) const;
+
+        //- Update result field based on interface functionality
+        virtual void updateInterfaceMatrix
+        (
+            solveScalarField& result,
+            const bool add,
+            const lduAddressing& lduAddr,
+            const label patchId,
+            const solveScalarField& psiInternal,
+            const scalarField& coeffs,
+            const direction cmpt,
+            const Pstream::commsTypes commsType
+        ) const;
+
+        //- Initialise neighbour matrix update
+        virtual void initInterfaceMatrixUpdate
+        (
+            Field<scalar>& result,
+            const bool add,
+            const lduAddressing& lduAddr,
+            const label patchId,
+            const Field<scalar>& psiInternal,
+            const scalarField& coeffs,
+            const Pstream::commsTypes commsType
+        ) const
+        {
+            NotImplemented;
+        }
+
+        //- Update result field based on interface functionality
+        virtual void updateInterfaceMatrix
+        (
+            Field<scalar>& result,
+            const bool add,
+            const lduAddressing& lduAddr,
+            const label patchId,
+            const Field<scalar>& psiInternal,
+            const scalarField& coeffs,
+            const Pstream::commsTypes commsType
+        ) const
+        {
+            NotImplemented;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "lduCalculatedProcessorField.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduPrimitiveProcessorInterface.C b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduPrimitiveProcessorInterface.C
index 2df7eeec30..48074e2e3e 100644
--- a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduPrimitiveProcessorInterface.C
+++ b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduPrimitiveProcessorInterface.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -38,6 +38,20 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::lduPrimitiveProcessorInterface::lduPrimitiveProcessorInterface
+(
+    const lduPrimitiveProcessorInterface& p
+)
+:
+    faceCells_(p.faceCells()),
+    myProcNo_(p.myProcNo()),
+    neighbProcNo_(p.neighbProcNo()),
+    forwardT_(p.forwardT()),
+    tag_(p.tag()),
+    comm_(p.comm())
+{}
+
+
 Foam::lduPrimitiveProcessorInterface::lduPrimitiveProcessorInterface
 (
     const labelUList& faceCells,
diff --git a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduPrimitiveProcessorInterface.H b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduPrimitiveProcessorInterface.H
index 47217dc6a2..f13a2b7f96 100644
--- a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduPrimitiveProcessorInterface.H
+++ b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduPrimitiveProcessorInterface.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -78,11 +78,6 @@ class lduPrimitiveProcessorInterface
 
     // Private Member Functions
 
-        //- No copy construct
-        lduPrimitiveProcessorInterface
-        (
-            const lduPrimitiveProcessorInterface&
-        ) = delete;
 
         //- No copy assignment
         void operator=(const lduPrimitiveProcessorInterface&) = delete;
@@ -107,6 +102,11 @@ public:
             const label comm = UPstream::worldComm
         );
 
+        //- Copy constructor
+        lduPrimitiveProcessorInterface
+        (
+            const lduPrimitiveProcessorInterface&
+        );
 
     //- Destructor
     virtual ~lduPrimitiveProcessorInterface() = default;
diff --git a/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C b/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
index 7042b98f62..cacc91b306 100644
--- a/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
+++ b/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
@@ -33,6 +33,7 @@ License
 #include "typeInfo.H"
 #include "addToRunTimeSelectionTable.H"
 #include "boundaryRadiationProperties.H"
+#include "lduCalculatedProcessorField.H"
 
 using namespace Foam::constant;
 
@@ -77,6 +78,8 @@ void Foam::radiation::viewFactor::initialise()
     DebugInFunction
         << "Total number of clusters : " << totalNCoarseFaces_ << endl;
 
+    useDirect_ = coeffs_.getOrDefault<bool>("useDirectSolver", true);
+
     map_.reset
     (
         new IOmapDistribute
@@ -93,94 +96,385 @@ void Foam::radiation::viewFactor::initialise()
         )
     );
 
-    scalarListIOList FmyProc
+    FmyProc_.reset
     (
-        IOobject
+        new scalarListIOList
         (
-            "F",
-            mesh_.facesInstance(),
-            mesh_,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE,
-            false
+            IOobject
+            (
+                "F",
+                mesh_.facesInstance(),
+                mesh_,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                false
+            )
         )
     );
 
-    labelListIOList globalFaceFaces
+    globalFaceFaces_.reset
     (
-        IOobject
+        new labelListIOList
         (
-            "globalFaceFaces",
-            mesh_.facesInstance(),
-            mesh_,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE,
-            false
+            IOobject
+            (
+                "globalFaceFaces",
+                mesh_.facesInstance(),
+                mesh_,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                false
+            )
         )
     );
 
-    List<labelListList> globalFaceFacesProc(Pstream::nProcs());
-    globalFaceFacesProc[Pstream::myProcNo()] = globalFaceFaces;
-    Pstream::gatherList(globalFaceFacesProc);
-
-    List<scalarListList> F(Pstream::nProcs());
-    F[Pstream::myProcNo()] = FmyProc;
-    Pstream::gatherList(F);
 
     globalIndex globalNumbering(nLocalCoarseFaces_);
 
-    if (Pstream::master())
+    // Size coarse qr
+    qrBandI_.setSize(nBands_);
+    forAll (qrBandI_, bandI)
     {
-        Fmatrix_.reset
-        (
-            new scalarSquareMatrix(totalNCoarseFaces_, Zero)
-        );
+        qrBandI_[bandI].setSize(nLocalCoarseFaces_, 0.0);
+    }
 
-        DebugInFunction
-            << "Insert elements in the matrix..." << endl;
+    if (!useDirect_)
+    {
+        DynamicList<label> dfaceJ;
 
-        for (const int procI : Pstream::allProcs())
+        // Per processor to owner (local)/neighbour (remote)
+        List<DynamicList<label>> procOwner(Pstream::nProcs());
+        List<DynamicList<label>> dynProcNeighbour(Pstream::nProcs());
+
+        forAll (globalFaceFaces_(), iFace)
         {
-            insertMatrixElements
-            (
-                globalNumbering,
-                procI,
-                globalFaceFacesProc[procI],
-                F[procI],
-                Fmatrix_()
-            );
-        }
-
-
-        if (coeffs_.get<bool>("smoothing"))
-        {
-            DebugInFunction << "Smoothing the matrix..." << endl;
-
-            for (label i=0; i<totalNCoarseFaces_; i++)
+            const labelList& visFaces = globalFaceFaces_()[iFace];
+            forAll (visFaces, j)
             {
-                scalar sumF = 0.0;
-                for (label j=0; j<totalNCoarseFaces_; j++)
-                {
-                    sumF += Fmatrix_()(i, j);
-                }
+                label gFacej = visFaces[j];
+                label proci = globalNumbering.whichProcID(gFacej);
+                label faceJ = globalNumbering.toLocal(proci, gFacej);
 
-                const scalar delta = sumF - 1.0;
-                for (label j=0; j<totalNCoarseFaces_; j++)
+                if (Pstream::myProcNo() == proci)
                 {
-                    Fmatrix_()(i, j) *= (1.0 - delta/(sumF + 0.001));
+                    edge e(iFace, faceJ);
+                    if (rays_.insert(e))
+                    {
+                        dfaceJ.append(j);
+                    }
+                }
+                else
+                {
+                    label gFaceI =
+                        globalNumbering.toGlobal(Pstream::myProcNo(), iFace);
+
+                    label proci = globalNumbering.whichProcID(gFacej);
+
+                    label facei =
+                        globalNumbering.toLocal(Pstream::myProcNo(), gFaceI);
+
+                    label remoteFacei = globalNumbering.toLocal(proci, gFacej);
+
+                    procOwner[proci].append(facei);
+                    dynProcNeighbour[proci].append(remoteFacei);
                 }
             }
         }
 
-        coeffs_.readEntry("constantEmissivity", constEmissivity_);
-        if (constEmissivity_)
+        mapRayToFmy_.transfer(dfaceJ);
+
+        labelList upper(rays_.size(), -1);
+        labelList lower(rays_.size(), -1);
+
+        const edgeList& raysLst = rays_.sortedToc();
+        label rayI = 0;
+        for (const auto& e : raysLst)
         {
-            CLU_.reset
+            label faceI = e.start();
+            label faceJ = e.end();
+            upper[rayI] = max(faceI, faceJ);
+            lower[rayI] = min(faceI, faceJ);
+            rayI++;
+        }
+
+        labelListList procNeighbour(dynProcNeighbour.size());
+        forAll(procNeighbour, i)
+        {
+            procNeighbour[i] = std::move(dynProcNeighbour[i]);
+        }
+        labelListList mySendCells;
+        Pstream::exchange<labelList, label>(procNeighbour, mySendCells);
+
+        label nbri = 0;
+        forAll(procOwner, proci)
+        {
+            if (procOwner[proci].size())
+            {
+                nbri++;
+            }
+            if (mySendCells[proci].size())
+            {
+                nbri++;
+            }
+        }
+
+        DebugInFunction<< "Number of procBound : " <<  nbri << endl;
+
+        PtrList<const lduPrimitiveProcessorInterface> primitiveInterfaces;
+
+        primitiveInterfaces.setSize(nbri);
+        internalCoeffs_.setSize(0);
+        boundaryCoeffs_.setSize(nbri);
+
+        nbri = 0;
+
+        procToInterface_.setSize(Pstream::nProcs(), -1);
+
+        forAll(procOwner, proci)
+        {
+            if (proci < Pstream::myProcNo() && procOwner[proci].size())
+            {
+                if (debug)
+                {
+                    Pout<< "Adding interface " << nbri
+                        << " to receive my " << procOwner[proci].size()
+                        << " from " << proci << endl;
+                }
+                procToInterface_[proci] = nbri;
+                primitiveInterfaces.set
+                (
+                    nbri++,
+                    new lduPrimitiveProcessorInterface
+                    (
+                        procOwner[proci],
+                        Pstream::myProcNo(),
+                        proci,
+                        tensorField(0),
+                        Pstream::msgType()+2
+                    )
+                );
+            }
+            else if (proci > Pstream::myProcNo() && mySendCells[proci].size())
+            {
+                if (debug)
+                {
+                    Pout<< "Adding interface " << nbri
+                        << " to send my " << mySendCells[proci].size()
+                        << " to " << proci << endl;
+                }
+                primitiveInterfaces.set
+                (
+                    nbri++,
+                    new lduPrimitiveProcessorInterface
+                    (
+                        mySendCells[proci],
+                        Pstream::myProcNo(),
+                        proci,
+                        tensorField(0),
+                        Pstream::msgType()+2
+                    )
+                );
+            }
+        }
+        forAll(procOwner, proci)
+        {
+            if (proci > Pstream::myProcNo() && procOwner[proci].size())
+            {
+                if (debug)
+                {
+                    Pout<< "Adding interface " << nbri
+                        << " to receive my " << procOwner[proci].size()
+                        << " from " << proci << endl;
+                }
+                procToInterface_[proci] = nbri;
+                primitiveInterfaces.set
+                (
+                    nbri++,
+                    new lduPrimitiveProcessorInterface
+                    (
+                        procOwner[proci],
+                        Pstream::myProcNo(),
+                        proci,
+                        tensorField(0),
+                        Pstream::msgType()+3
+                    )
+                );
+            }
+            else if (proci < Pstream::myProcNo() && mySendCells[proci].size())
+            {
+                if (debug)
+                {
+                    Pout<< "Adding interface " << nbri
+                        << " to send my " << mySendCells[proci].size()
+                        << " to " << proci << endl;
+                }
+                primitiveInterfaces.set
+                (
+                    nbri++,
+                    new lduPrimitiveProcessorInterface
+                    (
+                        mySendCells[proci],
+                        Pstream::myProcNo(),
+                        proci,
+                        tensorField(0),
+                        Pstream::msgType()+3
+                    )
+                );
+            }
+        }
+
+        forAll (boundaryCoeffs_, proci)
+        {
+            boundaryCoeffs_.set
+            (
+                proci,
+                new Field<scalar>
+                (
+                    primitiveInterfaces[proci].faceCells().size(),
+                    Zero
+                )
+            );
+        }
+
+        lduInterfacePtrsList allInterfaces;
+        allInterfaces.setSize(primitiveInterfaces.size());
+
+        forAll(primitiveInterfaces, i)
+        {
+            const lduPrimitiveProcessorInterface& pp = primitiveInterfaces[i];
+
+            allInterfaces.set(i, &pp);
+        }
+
+        const lduSchedule ps
+        (
+            lduPrimitiveMesh::nonBlockingSchedule
+                <lduPrimitiveProcessorInterface>(allInterfaces)
+        );
+
+        PtrList<const lduInterface> allInterfacesPtr(allInterfaces.size());
+        forAll (allInterfacesPtr, i)
+        {
+            const lduPrimitiveProcessorInterface& pp = primitiveInterfaces[i];
+
+            allInterfacesPtr.set
+            (
+                i,
+                new lduPrimitiveProcessorInterface(pp)
+            );
+        }
+
+        lduPtr_.reset
+        (
+            new lduPrimitiveMesh
+            (
+                rays_.size(),
+                lower,
+                upper,
+                allInterfacesPtr,
+                ps,
+                UPstream::worldComm
+            )
+        );
+
+        // Set size for local lduMatrix
+        matrixPtr_.reset(new lduMatrix(lduPtr_()));
+
+        scalarListList& myF = FmyProc_();
+
+        if (coeffs_.get<bool>("smoothing"))
+        {
+            scalar maxDelta = 0;
+            scalar totalDelta = 0;
+            forAll (myF, i)
+            {
+                scalar sumF = 0.0;
+                scalarList& myFij = myF[i];
+                forAll (myFij, j)
+                {
+                    sumF += myFij[j];
+                }
+                const scalar delta = sumF - 1.0;
+                forAll (myFij, j)
+                {
+                    myFij[j] *= (1.0 - delta/(sumF + 0.001));
+                }
+                totalDelta += delta;
+                if (delta > maxDelta)
+                {
+                    maxDelta = delta;
+                }
+            }
+            totalDelta /= myF.size();
+            reduce(totalDelta, sumOp<scalar>());
+            reduce(maxDelta, maxOp<scalar>());
+            Info << "Smoothng average delta : " << totalDelta << endl;
+            Info << "Smoothng maximum delta : " << maxDelta << nl << endl;
+        }
+    }
+
+    if (useDirect_)
+    {
+        List<labelListList> globalFaceFacesProc(Pstream::nProcs());
+        globalFaceFacesProc[Pstream::myProcNo()] = globalFaceFaces_();
+        Pstream::gatherList(globalFaceFacesProc);
+
+        List<scalarListList> F(Pstream::nProcs());
+        F[Pstream::myProcNo()] = FmyProc_();
+        Pstream::gatherList(F);
+
+        if (Pstream::master())
+        {
+            Fmatrix_.reset
             (
                 new scalarSquareMatrix(totalNCoarseFaces_, Zero)
             );
 
-            pivotIndices_.setSize(CLU_().m());
+            DebugInFunction
+                << "Insert elements in the matrix..." << endl;
+
+            for (const int procI : Pstream::allProcs())
+            {
+                insertMatrixElements
+                (
+                    globalNumbering,
+                    procI,
+                    globalFaceFacesProc[procI],
+                    F[procI],
+                    Fmatrix_()
+                );
+            }
+
+            if (coeffs_.get<bool>("smoothing"))
+            {
+                DebugInFunction << "Smoothing the matrix..." << endl;
+
+                for (label i=0; i<totalNCoarseFaces_; i++)
+                {
+                    scalar sumF = 0.0;
+                    for (label j=0; j<totalNCoarseFaces_; j++)
+                    {
+                        sumF += Fmatrix_()(i, j);
+                    }
+
+                    const scalar delta = sumF - 1.0;
+                    for (label j=0; j<totalNCoarseFaces_; j++)
+                    {
+                        Fmatrix_()(i, j) *= (1.0 - delta/(sumF + 0.001));
+                    }
+                }
+            }
+
+            coeffs_.readEntry("constantEmissivity", constEmissivity_);
+            if (constEmissivity_)
+            {
+                CLU_.reset
+                (
+                    new scalarSquareMatrix(totalNCoarseFaces_, Zero)
+                );
+
+                pivotIndices_.setSize(CLU_().m());
+            }
         }
     }
 
@@ -258,7 +552,8 @@ Foam::radiation::viewFactor::viewFactor(const volScalarField& T)
     pivotIndices_(0),
     useSolarLoad_(false),
     solarLoad_(),
-    nBands_(coeffs_.getOrDefault<label>("nBands", 1))
+    nBands_(coeffs_.getOrDefault<label>("nBands", 1)),
+    FmyProc_()
 {
     initialise();
 }
@@ -320,7 +615,8 @@ Foam::radiation::viewFactor::viewFactor
     pivotIndices_(0),
     useSolarLoad_(false),
     solarLoad_(),
-    nBands_(coeffs_.getOrDefault<label>("nBands", 1))
+    nBands_(coeffs_.getOrDefault<label>("nBands", 1)),
+    FmyProc_()
 {
     initialise();
 }
@@ -372,8 +668,13 @@ void Foam::radiation::viewFactor::calculate()
         solarLoad_->calculate();
     }
 
-     // Net radiation
-    scalarField q(totalNCoarseFaces_, 0.0);
+    // Global net radiation
+    scalarField qNet(totalNCoarseFaces_, 0.0);
+
+    // Local net radiation
+    scalarField qTotalCoarse(nLocalCoarseFaces_, 0.0);
+
+    // Referen to fvMesh qr field
     volScalarField::Boundary& qrBf = qr_.boundaryFieldRef();
 
     globalIndex globalNumbering(nLocalCoarseFaces_);
@@ -383,6 +684,9 @@ void Foam::radiation::viewFactor::calculate()
 
     for (label bandI = 0; bandI < nBands_; bandI++)
     {
+        // Global bandI radiation
+        scalarField qBandI(totalNCoarseFaces_, 0.0);
+
         scalarField compactCoarseT4(map_->constructSize(), 0.0);
         scalarField compactCoarseE(map_->constructSize(), 0.0);
         scalarField compactCoarseHo(map_->constructSize(), 0.0);
@@ -467,6 +771,7 @@ void Foam::radiation::viewFactor::calculate()
         SubList<scalar>(compactCoarseHo, nLocalCoarseFaces_) =
             localCoarseHoave;
 
+
         // Distribute data
         map_->distribute(compactCoarseT4);
         map_->distribute(compactCoarseE);
@@ -487,123 +792,260 @@ void Foam::radiation::viewFactor::calculate()
 
         map_->distribute(compactGlobalIds);
 
-        // Create global size vectors
-        scalarField T4(totalNCoarseFaces_, 0.0);
-        scalarField E(totalNCoarseFaces_, 0.0);
-        scalarField qrExt(totalNCoarseFaces_, 0.0);
-
-        // Fill lists from compact to global indexes.
-        forAll(compactCoarseT4, i)
+        if (!useDirect_)
         {
-            T4[compactGlobalIds[i]] = compactCoarseT4[i];
-            E[compactGlobalIds[i]] = compactCoarseE[i];
-            qrExt[compactGlobalIds[i]] = compactCoarseHo[i];
-        }
+            const labelList globalToCompact
+            (
+                invert(totalNCoarseFaces_, compactGlobalIds)
+            );
 
-        Pstream::listCombineAllGather(T4, maxEqOp<scalar>());
-        Pstream::listCombineAllGather(E, maxEqOp<scalar>());
-        Pstream::listCombineAllGather(qrExt, maxEqOp<scalar>());
+            scalarField& diag = matrixPtr_->diag(localCoarseEave.size());
+            scalarField& upper = matrixPtr_->upper(rays_.size());
+            scalarField& lower = matrixPtr_->lower(rays_.size());
 
-        if (Pstream::master())
-        {
-            // Variable emissivity
-            if (!constEmissivity_)
+            scalarField source(nLocalCoarseFaces_, 0);
+
+
+            // Local diag and source
+            forAll(source, i)
             {
-                scalarSquareMatrix C(totalNCoarseFaces_, 0.0);
+                const scalar sigmaT4 =
+                    physicoChemical::sigma.value()*localCoarseT4ave[i];
 
-                for (label i=0; i<totalNCoarseFaces_; i++)
+                diag[i] = 1/localCoarseEave[i];
+
+                source[i] += -sigmaT4 + localCoarseHoave[i];
+            }
+
+            // Local matrix coefficients
+            if (!constEmissivity_ || iterCounter_ == 0)
+            {
+                const edgeList& raysLst = rays_.sortedToc();
+
+                label rayI = 0;
+                for (const auto& e : raysLst)
                 {
-                    for (label j=0; j<totalNCoarseFaces_; j++)
+                    label facelJ = e.end();
+                    label faceI = e.start();
+
+                    label faceJ = mapRayToFmy_[rayI];
+
+                    lower[rayI] =
+                        (1-1/localCoarseEave[faceI])*FmyProc_()[faceI][faceJ];
+
+                    upper[rayI] =
+                        (1-1/localCoarseEave[facelJ])*FmyProc_()[faceI][faceJ];
+
+                    rayI++;
+                }
+                iterCounter_++;
+            }
+
+            // Extra local contribution to the source and extra matrix coefs
+            // from other procs (boundaryCoeffs)
+            label nInterfaces = lduPtr_().interfaces().size();
+            labelList boundCoeffI(nInterfaces, Zero);
+            forAll (globalFaceFaces_(), iFace)
+            {
+                const labelList& visFaces = globalFaceFaces_()[iFace];
+                forAll (visFaces, jFace)
+                {
+                    label gFacej = visFaces[jFace];
+                    label proci = globalNumbering.whichProcID(gFacej);
+                    if (Pstream::myProcNo() == proci)
                     {
-                        const scalar invEj = 1.0/E[j];
+                        label lFacej =
+                            globalNumbering.toLocal
+                            (
+                                Pstream::myProcNo(),
+                                gFacej
+                            );
+
                         const scalar sigmaT4 =
-                            physicoChemical::sigma.value()*T4[j];
+                            physicoChemical::sigma.value()
+                            *localCoarseT4ave[lFacej];
 
-                        if (i==j)
-                        {
-                            C(i, j) = invEj - (invEj - 1.0)*Fmatrix_()(i, j);
-                            q[i] +=
-                                (Fmatrix_()(i, j) - 1.0)*sigmaT4 + qrExt[j];
-                        }
-                        else
-                        {
-                            C(i, j) = (1.0 - invEj)*Fmatrix_()(i, j);
-                            q[i] += Fmatrix_()(i, j)*sigmaT4;
-                        }
+                        source[iFace] += FmyProc_()[iFace][jFace]*sigmaT4;
+                    }
+                    else
+                    {
+                        label compactFaceJ = globalToCompact[gFacej];
+                        const scalar sigmaT4 =
+                            physicoChemical::sigma.value()
+                            *compactCoarseT4[compactFaceJ];
 
+                        source[iFace] += FmyProc_()[iFace][jFace]*sigmaT4;
+
+                        label interfaceI = procToInterface_[proci];
+
+                        boundaryCoeffs_
+                            [interfaceI][boundCoeffI[interfaceI]++] =
+                                -(1-1/compactCoarseE[compactFaceJ])
+                                *FmyProc_()[iFace][jFace];
                     }
                 }
-
-                Info<< "Solving view factor equations for band :"
-                    << bandI << endl;
-
-                // Negative coming into the fluid
-                LUsolve(C, q);
             }
-            else //Constant emissivity
+
+            PtrList<const lduInterfaceField> interfaces(nInterfaces);
+            for(label i = 0; i < interfaces.size(); i++)
             {
-                // Initial iter calculates CLU and caches it
-                if (iterCounter_ == 0)
+                interfaces.set
+                (
+                    i,
+                    new lduCalculatedProcessorField<scalar>
+                    (
+                        lduPtr_().interfaces()[i],
+                        qrBandI_[bandI]
+                    )
+                );
+            }
+
+            const dictionary& solverControls =
+                qr_.mesh().solverDict
+                (
+                    qr_.select
+                    (
+                        qr_.mesh().data::template getOrDefault<bool>
+                        ("finalIteration", false)
+                    )
+                );
+
+            // Solver call
+            solverPerformance solverPerf = lduMatrix::solver::New
+            (
+                "qr",
+                matrixPtr_(),
+                boundaryCoeffs_,
+                internalCoeffs_,
+                interfaces,
+                solverControls
+            )->solve(qrBandI_[bandI], source);
+
+            solverPerf.print(Info.masterStream(qr_.mesh().comm()));
+
+            qTotalCoarse += qrBandI_[bandI];
+        }
+
+        if (useDirect_)
+        {
+            // Create global size vectors
+            scalarField T4(totalNCoarseFaces_, 0.0);
+            scalarField E(totalNCoarseFaces_, 0.0);
+            scalarField qrExt(totalNCoarseFaces_, 0.0);
+
+            // Fill lists from compact to global indexes.
+            forAll(compactCoarseT4, i)
+            {
+                T4[compactGlobalIds[i]] = compactCoarseT4[i];
+                E[compactGlobalIds[i]] = compactCoarseE[i];
+                qrExt[compactGlobalIds[i]] = compactCoarseHo[i];
+            }
+
+            Pstream::listCombineGather(T4, maxEqOp<scalar>());
+            Pstream::listCombineGather(E, maxEqOp<scalar>());
+            Pstream::listCombineGather(qrExt, maxEqOp<scalar>());
+
+            Pstream::listCombineScatter(T4);
+            Pstream::listCombineScatter(E);
+            Pstream::listCombineScatter(qrExt);
+
+            if (Pstream::master())
+            {
+                // Variable emissivity
+                if (!constEmissivity_)
                 {
+                    scalarSquareMatrix C(totalNCoarseFaces_, 0.0);
+
                     for (label i=0; i<totalNCoarseFaces_; i++)
                     {
                         for (label j=0; j<totalNCoarseFaces_; j++)
                         {
                             const scalar invEj = 1.0/E[j];
+                            const scalar sigmaT4 =
+                                physicoChemical::sigma.value()*T4[j];
+
                             if (i==j)
                             {
-                                CLU_()(i, j) =
-                                    invEj-(invEj-1.0)*Fmatrix_()(i, j);
+                                C(i, j) = invEj;
+                                qBandI[i] += -sigmaT4 +  qrExt[j];
                             }
                             else
                             {
-                                CLU_()(i, j) = (1.0 - invEj)*Fmatrix_()(i, j);
+                                C(i, j) = (1.0 - invEj)*Fmatrix_()(i, j);
+                                qBandI[i] += Fmatrix_()(i, j)*sigmaT4;
                             }
                         }
                     }
 
-                    DebugInFunction << "\nDecomposing C matrix..." << endl;
+                    Info<< "Solving view factor equations for band :"
+                        << bandI << endl;
 
-                    LUDecompose(CLU_(), pivotIndices_);
+                    // Negative coming into the fluid
+                    LUsolve(C, qBandI);
                 }
-
-                for (label i=0; i<totalNCoarseFaces_; i++)
+                else //Constant emissivity
                 {
-                    for (label j=0; j<totalNCoarseFaces_; j++)
+                    // Initial iter calculates CLU and caches it
+                    if (iterCounter_ == 0)
                     {
-                        const scalar sigmaT4 =
-                            constant::physicoChemical::sigma.value()*T4[j];
-
-                        if (i==j)
+                        for (label i=0; i<totalNCoarseFaces_; i++)
                         {
-                            q[i] +=
-                            (Fmatrix_()(i, j) - 1.0)*sigmaT4  + qrExt[j];
+                            for (label j=0; j<totalNCoarseFaces_; j++)
+                            {
+                                const scalar invEj = 1.0/E[j];
+                                if (i==j)
+                                {
+                                    CLU_()(i, j) = invEj;
+                                }
+                                else
+                                {
+                                    CLU_()(i, j) =
+                                        (1.0-invEj)*Fmatrix_()(i, j);
+                                }
+                            }
                         }
-                        else
+
+                        DebugInFunction << "\nDecomposing C matrix..." << endl;
+
+                        LUDecompose(CLU_(), pivotIndices_);
+                    }
+
+                    for (label i=0; i<totalNCoarseFaces_; i++)
+                    {
+                        for (label j=0; j<totalNCoarseFaces_; j++)
                         {
-                            q[i] += Fmatrix_()(i, j)*sigmaT4;
+                            const scalar sigmaT4 =
+                                constant::physicoChemical::sigma.value()*T4[j];
+
+                            if (i==j)
+                            {
+                                qBandI[i] += -sigmaT4 + qrExt[j];
+                            }
+                            else
+                            {
+                                qBandI[i] += Fmatrix_()(i, j)*sigmaT4;
+                            }
                         }
                     }
+
+                    Info<< "Solving view factor equations for band : "
+                        << bandI  << endl;
+
+                    LUBacksubstitute(CLU_(), pivotIndices_, qBandI);
+                    iterCounter_ ++;
                 }
-
-
-                Info<< "Solving view factor equations for band : "
-                    << bandI  << endl;
-
-
-                LUBacksubstitute(CLU_(), pivotIndices_, q);
-                iterCounter_ ++;
             }
-        }
+            // Broadcast qBandI and fill qr
+            Pstream::broadcast(qBandI);
 
+            qNet += qBandI;
+        }
     }
-    // Broadcast q and fill qr
-    Pstream::broadcast(q);
 
     label globCoarseId = 0;
     for (const label patchID : selectedPatches_)
     {
-        const polyPatch& pp = mesh_.boundaryMesh()[patchID];
+        const polyPatch& pp = coarseMesh_.boundaryMesh()[patchID];
 
         if (pp.size() > 0)
         {
@@ -617,22 +1059,31 @@ void Foam::radiation::viewFactor::calculate()
             const labelList& coarsePatchFace =
                 coarseMesh_.patchFaceMap()[patchID];
 
-            /// scalar heatFlux = 0.0;
+            //scalar heatFlux = 0.0;
             forAll(coarseToFine, coarseI)
             {
-                label globalCoarse =
-                    globalNumbering.toGlobal
-                    (Pstream::myProcNo(), globCoarseId);
+                label globalCoarse = globalNumbering.toGlobal
+                (
+                    Pstream::myProcNo(),
+                    globCoarseId
+                );
 
                 const label coarseFaceID = coarsePatchFace[coarseI];
                 const labelList& fineFaces = coarseToFine[coarseFaceID];
 
                 for (const label facei : fineFaces)
                 {
-                    qrp[facei] = q[globalCoarse];
-                    /// heatFlux += qrp[facei]*sf[facei];
+                    if (useDirect_)
+                    {
+                        qrp[facei] = qNet[globalCoarse];
+                    }
+                    else
+                    {
+                        qrp[facei] = qTotalCoarse[globCoarseId];
+                    }
+                    //heatFlux += qrp[facei]*sf[facei];
                 }
-                globCoarseId ++;
+                globCoarseId++;
             }
         }
     }
@@ -645,8 +1096,7 @@ void Foam::radiation::viewFactor::calculate()
             const scalarField& magSf = mesh_.magSf().boundaryField()[patchID];
             const scalar heatFlux = gSum(qrp*magSf);
 
-            InfoInFunction
-                << "Total heat transfer rate at patch: "
+            Info<< "Total heat transfer rate at patch: "
                 << patchID << " "
                 << heatFlux << endl;
         }
diff --git a/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.H b/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.H
index 7592cfdf50..662c6586be 100644
--- a/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.H
+++ b/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.H
@@ -60,6 +60,10 @@ SourceFiles
 #include "IOmapDistribute.H"
 #include "solarLoad.H"
 
+#include "lduPrimitiveMesh.H"
+#include "lduPrimitiveProcessorInterface.H"
+#include "lduMatrix.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -108,12 +112,18 @@ protected:
         //- Net radiative heat flux [W/m2]
         volScalarField qr_;
 
+        //- Coarse radiative heat flux
+        List<scalarField> qrBandI_;
+
         //- View factor matrix
         autoPtr<scalarSquareMatrix> Fmatrix_;
 
         //- Inverse of C matrix
         autoPtr<scalarSquareMatrix> CLU_;
 
+        //- Visible global faces
+        autoPtr<labelListIOList> globalFaceFaces_;
+
         //- Selected patches
         labelList selectedPatches_;
 
@@ -141,6 +151,35 @@ protected:
         //-Number of bands
         label nBands_;
 
+        //- Primitive addressing for lduMatrix
+        autoPtr<lduPrimitiveMesh> lduPtr_;
+
+        //- Matrix formed from view factors
+        autoPtr<lduMatrix> matrixPtr_;
+
+        //- Boundary scalar field containing pseudo-matrix coeffs
+        //- for internal cells
+        FieldField<Field, scalar> internalCoeffs_;
+
+        //- Boundary scalar field containing pseudo-matrix coeffs
+        //- for boundary cells
+        FieldField<Field, scalar> boundaryCoeffs_;
+
+        //- Rays on local proc
+        edgeHashSet rays_;
+
+        //- Map local-ray to j-column for F
+        List<label> mapRayToFmy_;
+
+        //- Local view factors
+        autoPtr<scalarListIOList> FmyProc_;
+
+        //- Map from proc to interafce
+        labelList procToInterface_;
+
+        //- Use direct or iterative solver
+        bool useDirect_;
+
 
     // Private Member Functions