From 913c45affff2144f511adce6474548ebee0fdb45 Mon Sep 17 00:00:00 2001
From: Andrew Heather <>
Date: Fri, 9 Sep 2022 13:12:15 +0100
Subject: [PATCH] ENH: faMatrix - added dot operator functions

---
 src/finiteArea/faMatrices/faMatrix/faMatrix.C | 121 ++++++++++++++++++
 src/finiteArea/faMatrices/faMatrix/faMatrix.H |  72 +++++++++++
 2 files changed, 193 insertions(+)

diff --git a/src/finiteArea/faMatrices/faMatrix/faMatrix.C b/src/finiteArea/faMatrices/faMatrix/faMatrix.C
index 2307a1bb79..e5a1eed1a5 100644
--- a/src/finiteArea/faMatrices/faMatrix/faMatrix.C
+++ b/src/finiteArea/faMatrices/faMatrix/faMatrix.C
@@ -2074,6 +2074,127 @@ Foam::tmp<Foam::faMatrix<Type>> Foam::operator*
 }
 
 
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::faPatchField, Foam::areaMesh>>
+Foam::operator&
+(
+    const faMatrix<Type>& M,
+    const DimensionedField<Type, areaMesh>& psi
+)
+{
+    auto tMphi = tmp<GeometricField<Type, faPatchField, areaMesh>>::New
+    (
+        IOobject
+        (
+            "M&" + psi.name(),
+            psi.instance(),
+            psi.mesh().mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        psi.mesh(),
+        M.dimensions()/dimArea
+    );
+    auto& Mphi = tMphi.ref();
+
+    // Loop over field components
+    if (M.hasDiag())
+    {
+        for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
+        {
+            scalarField psiCmpt(psi.field().component(cmpt));
+            scalarField boundaryDiagCmpt(M.diag());
+            M.addBoundaryDiag(boundaryDiagCmpt, cmpt);
+            Mphi.primitiveFieldRef().replace(cmpt, -boundaryDiagCmpt*psiCmpt);
+        }
+    }
+    else
+    {
+        Mphi.primitiveFieldRef() = Zero;
+    }
+
+    Mphi.primitiveFieldRef() += M.lduMatrix::H(psi.field()) + M.source();
+    M.addBoundarySource(Mphi.primitiveFieldRef());
+
+    Mphi.primitiveFieldRef() /= -psi.mesh().S();
+    Mphi.correctBoundaryConditions();
+
+    return tMphi;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::faPatchField, Foam::areaMesh>>
+Foam::operator&
+(
+    const faMatrix<Type>& M,
+    const tmp<DimensionedField<Type, areaMesh>>& tpsi
+)
+{
+    tmp<GeometricField<Type, faPatchField, areaMesh>> tMpsi = M & tpsi();
+    tpsi.clear();
+    return tMpsi;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::faPatchField, Foam::areaMesh>>
+Foam::operator&
+(
+    const faMatrix<Type>& M,
+    const tmp<GeometricField<Type, faPatchField, areaMesh>>& tpsi
+)
+{
+    tmp<GeometricField<Type, faPatchField, areaMesh>> tMpsi = M & tpsi();
+    tpsi.clear();
+    return tMpsi;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::faPatchField, Foam::areaMesh>>
+Foam::operator&
+(
+    const tmp<faMatrix<Type>>& tM,
+    const DimensionedField<Type, areaMesh>& psi
+)
+{
+    tmp<GeometricField<Type, faPatchField, areaMesh>> tMpsi = tM() & psi;
+    tM.clear();
+    return tMpsi;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::faPatchField, Foam::areaMesh>>
+Foam::operator&
+(
+    const tmp<faMatrix<Type>>& tM,
+    const tmp<DimensionedField<Type, areaMesh>>& tpsi
+)
+{
+    tmp<GeometricField<Type, faPatchField, areaMesh>> tMpsi = tM() & tpsi();
+    tM.clear();
+    tpsi.clear();
+    return tMpsi;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::faPatchField, Foam::areaMesh>>
+Foam::operator&
+(
+    const tmp<faMatrix<Type>>& tM,
+    const tmp<GeometricField<Type, faPatchField, areaMesh>>& tpsi
+)
+{
+    tmp<GeometricField<Type, faPatchField, areaMesh>> tMpsi = tM() & tpsi();
+    tM.clear();
+    tpsi.clear();
+    return tMpsi;
+}
+
+
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
 template<class Type>
diff --git a/src/finiteArea/faMatrices/faMatrix/faMatrix.H b/src/finiteArea/faMatrices/faMatrix/faMatrix.H
index b0fc26c2e1..e5b6dfae97 100644
--- a/src/finiteArea/faMatrices/faMatrix/faMatrix.H
+++ b/src/finiteArea/faMatrices/faMatrix/faMatrix.H
@@ -64,6 +64,47 @@ template<class T> class UIndirectList;
 template<class Type>
 Ostream& operator<<(Ostream&, const faMatrix<Type>&);
 
+template<class Type>
+tmp<GeometricField<Type, faPatchField, areaMesh>> operator&
+(
+    const faMatrix<Type>&,
+    const DimensionedField<Type, areaMesh>&
+);
+
+template<class Type>
+tmp<GeometricField<Type, faPatchField, areaMesh>> operator&
+(
+    const faMatrix<Type>&,
+    const tmp<DimensionedField<Type, areaMesh>>&
+);
+
+template<class Type>
+tmp<GeometricField<Type, faPatchField, areaMesh>> operator&
+(
+    const faMatrix<Type>&,
+    const tmp<GeometricField<Type, faPatchField, areaMesh>>&
+);
+
+template<class Type>
+tmp<GeometricField<Type, faPatchField, areaMesh>> operator&
+(
+    const tmp<faMatrix<Type>>&,
+    const DimensionedField<Type, areaMesh>&
+);
+
+template<class Type>
+tmp<GeometricField<Type, faPatchField, areaMesh>> operator&
+(
+    const tmp<faMatrix<Type>>&,
+    const tmp<DimensionedField<Type, areaMesh>>&
+);
+
+template<class Type>
+tmp<GeometricField<Type, faPatchField, areaMesh>> operator&
+(
+    const tmp<faMatrix<Type>>&,
+    const tmp<GeometricField<Type, faPatchField, areaMesh>>&
+);
 
 /*---------------------------------------------------------------------------*\
                           Class faMatrix Declaration
@@ -469,6 +510,37 @@ public:
         void operator*=(const dimensioned<scalar>&);
 
 
+    // Friend Operators
+
+        friend tmp<GeometricField<Type, faPatchField, areaMesh>>
+        operator& <Type>
+        (
+            const faMatrix<Type>&,
+            const DimensionedField<Type, areaMesh>&
+        );
+
+        friend tmp<GeometricField<Type, faPatchField, areaMesh>>
+        operator& <Type>
+        (
+            const faMatrix<Type>&,
+            const tmp<GeometricField<Type, faPatchField, areaMesh>>&
+        );
+
+        friend tmp<GeometricField<Type, faPatchField, areaMesh>>
+        operator& <Type>
+        (
+            const tmp<faMatrix<Type>>&,
+            const DimensionedField<Type, areaMesh>&
+        );
+
+        friend tmp<GeometricField<Type, faPatchField, areaMesh>>
+        operator& <Type>
+        (
+            const tmp<faMatrix<Type>>&,
+            const tmp<GeometricField<Type, faPatchField, areaMesh>>&
+        );
+
+
     // Ostream Operator
 
         friend Ostream& operator<< <Type>