STYLE: ListExpression: copyright, indenting

Feature grad update

See merge request Development/openfoam!764

applications/test/expressionTemplates/Make/files

@@ -0,0 +1,5 @@
+Test-expressionTemplates-volFields.C
+/* Test-expressionTemplates-pointFields.C */
+/* Test-expressionTemplates-Fields.C */
+
+EXE = $(FOAM_USER_APPBIN)/Test-expressionTemplates

applications/test/expressionTemplates/Make/options

@@ -0,0 +1,11 @@
+EXE_INC = \
+    -O0 \
+    -I$(FOAM_SRC)/meshTools/lnInclude \
+    -I$(FOAM_SRC)/finiteVolume/lnInclude \
+
+EXE_LIBS = \
+/*    -rpath /Volumes/case_sensitive/OpenFOAM/work/feature-et/platforms */ \
+/*    -rpath /Volumes/case_sensitive/OpenFOAM/work/feature-et/platforms/darwin64GccDPInt32Opt/lib */ \
+/*    -rpath /Volumes/case_sensitive/OpenFOAM/work/feature-et/platforms/darwin64GccDPInt32Opt/lib/$(FOAM_MPI)  */\
+    -L/Volumes/case_sensitive/OpenFOAM/work/feature-et/platforms/darwin64GccDPInt32Opt/lib \
+    -lfiniteVolume

applications/test/expressionTemplates/Test-expressionTemplates-pointFields.C

@@ -0,0 +1,126 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2025 M. Janssens
+-------------------------------------------------------------------------------
+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/>.
+
+Application
+    Test-expressionTemplates
+
+\*---------------------------------------------------------------------------*/
+
+#include "Time.H"
+#include "argList.H"
+#include "polyMesh.H"
+#include "pointMesh.H"
+#include "pointFields.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+
+    #include "createTime.H"
+    #include "createPolyMesh.H"
+
+    {
+        scalarField vals0({1.0, 2.0, 3.0});
+        const Expression::ListConstRefWrap<scalar> wvals0(vals0);
+
+        scalarField vals1;
+        Expression::ListRefWrap<scalar> wvals1(vals1.size(), vals1);
+        wvals1 = wvals0;
+        return 0;
+    }
+
+
+
+
+    const pointMesh& pMesh = pointMesh::New(mesh);
+
+    // Field, dimensionedScalar as List expression
+    {
+        scalarField vals({1.0, 2.0, 3.0});
+        dimensionedScalar d(dimless, 4.0);
+
+        scalarField result(d.expr(vals.size())*vals.expr());
+
+        Pout<< "result:" << result << endl;
+    }
+
+    Info<< "Reading field p\n" << endl;
+    pointScalarField p
+    (
+        IOobject
+        (
+            "pointDisplacement",
+            runTime.timeName(),
+            pMesh.thisDb(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        pMesh
+    );
+
+    pointScalarField result
+    (
+        IOobject
+        (
+            "result",
+            runTime.timeName(),
+            pMesh.thisDb(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            IOobject::NO_REGISTER
+        ),
+        pMesh,
+        dimensionedScalar("zero", Foam::sqr(p().dimensions()), 0)
+    );
+
+    // Construct value + dimensions with correct sizing
+    const auto oneField(dimensionedScalar(dimless, 1.0).expr(p));
+
+    // Make expression
+    auto expression = oneField * p.expr() + p.expr();
+
+    // Combine expressions
+    auto newExpression = sqr(expression);
+
+    // Assign values
+    result = newExpression;
+
+    DebugVar(result);
+
+
+    // Construct from expression
+    pointScalarField result2("result2", pMesh, sqrt(p.expr()));
+
+    DebugVar(result2);
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/test/expressionTemplates/Test-expressionTemplates-volFields.C

@@ -0,0 +1,451 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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/>.
+
+Application
+    volPointInterpolationTest
+
+\*---------------------------------------------------------------------------*/
+
+#include "Time.H"
+#include "argList.H"
+#include "fvMesh.H"
+#include "ListExpression.H"
+#include "GeometricFieldExpression.H"
+#include "fvCFD.H"
+#include "fvMatrixExpression.H"
+#include <ratio>
+#include <chrono>
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+tmp<volScalarField> someFunction(const volScalarField& fld)
+{
+    return fld*1.0001;
+}
+
+
+template<class Type>
+void fusedGaussFvmLaplacian
+(
+    fvMatrix<Type>& fvm,
+    const surfaceInterpolationScheme<scalar>& interpGammaScheme,
+    const fv::snGradScheme<Type>& snGradScheme,
+    const GeometricField<scalar, fvPatchField, volMesh>& gamma,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    // Replacement for gaussLaplacianScheme::fvmLaplacian with scalar gamma
+    typedef GeometricField<Type, fvsPatchField, surfaceMesh> surfaceType;
+
+    const auto& mesh = vf.mesh();
+
+    // Expression for weights
+    const auto weights = interpGammaScheme.weights(gamma).expr();
+
+    // Expression for gamma_face * magSf
+    const auto gammaMagSf =
+        Expression::interpolate(gamma.expr(), weights, mesh)
+      * mesh.magSf().expr();
+
+    // Expression for deltaCoeffs
+    const auto deltaCoeffs = snGradScheme.deltaCoeffs(vf).expr();
+
+    // Construct matrix
+    Expression::fvmLaplacianUncorrected(fvm, gammaMagSf, deltaCoeffs);
+
+    if (snGradScheme.corrected())
+    {
+        // Wrap correction
+        const auto corr(snGradScheme.correction(vf).expr());
+        const auto V = mesh.V().expr();
+
+        if (mesh.fluxRequired(vf.name()))
+        {
+            fvm.faceFluxCorrectionPtr() = std::make_unique<surfaceType>
+            (
+                IOobject
+                (
+                    "faceFluxCorr",
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    IOobject::NO_REGISTER
+                ),
+                mesh,
+                gamma.dimensions()
+               *mesh.magSf().dimensions()
+               *corr.data().dimensions()
+            );
+            auto& faceFluxCorr = *fvm.faceFluxCorrectionPtr();
+            faceFluxCorr = gammaMagSf*corr;
+
+            fvm.source() =
+                fvm.source().expr()
+              - (
+                    V * fvc::div
+                    (
+                        faceFluxCorr
+                    )().primitiveField().expr()
+                );
+        }
+        else
+        {
+            // Temporary field
+            surfaceType faceFluxCorr
+            (
+                IOobject
+                (
+                    "faceFluxCorr",
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    IOobject::NO_REGISTER
+                ),
+                mesh,
+                gamma.dimensions()
+               *mesh.magSf().dimensions()
+               *corr.data().dimensions()
+            );
+            faceFluxCorr = gammaMagSf*corr;
+
+            fvm.source() =
+                fvm.source().expr()
+              - (
+                    V * fvc::div
+                    (
+                        faceFluxCorr
+                    )().primitiveField().expr()
+                );
+        }
+    }
+}
+
+
+using namespace std::chrono;
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh.thisDb(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+    
+    {
+        volScalarField p2("p2", p);
+        DebugVar(p2.boundaryFieldRef());
+    
+        for (auto& pf : p.boundaryFieldRef())
+        {
+            scalarField newVals(pf.size());
+            forAll(pf, i)
+            {
+                newVals[i] = scalar(i);
+            }
+            pf == newVals;
+        }
+
+        //std::vector<const scalarField*> ptrs;
+        // List<const scalarField*> ptrs;
+        // for (const auto& pp : p.boundaryField())
+        // {
+        //     ptrs.push_back(&pp);
+        // }
+        DebugVar(sizeof(long));
+        DebugVar(p.boundaryField());
+        Expression::ListsConstRefWrap<scalarField> expr(p.boundaryField());
+
+        const auto twoA = expr + expr;
+        Expression::ListsRefWrap<scalarField> expr2(p2.boundaryFieldRef());
+        Pout<< "**before assignment twoA:" << twoA.size()
+            << " expr2:" << expr2.size() << endl;
+        expr2 = twoA;
+        Pout<< "**after assignment twoA:" << twoA.size()
+            << " expr2:" << expr2.size() << endl;
+        DebugVar(p2.boundaryField());
+        // forAll(expr, i)
+        // {
+        //     Pout<< "i:" << i
+        //         //<< " expr:" << expr[i]
+        //         //<< " twoA:" << twoA[i]
+        //         << " expr2:" << expr2[i]
+        //         << endl;
+        // }
+        return 0;
+    }
+
+    {
+        //DebugVar(linearInterpolate(p));
+        //auto tweights = linear<scalar>(mesh).weights(p);
+        //DebugVar(tweights);
+
+        surfaceScalarField result
+        (
+            IOobject
+            (
+                "result",
+                runTime.timeName(),
+                mesh.thisDb(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                IOobject::NO_REGISTER
+            ),
+            mesh,
+            dimensionedScalar(p.dimensions(), 0)
+        );
+        //result = Expression::interpolate
+        //(
+        //    p.expr(),
+        //    tweights().expr(),
+        //    mesh
+        //);
+
+        result = Expression::linearInterpolate(p.expr(), mesh);
+
+
+        DebugVar(result);
+
+        return 0;
+    }
+
+
+
+    volScalarField p2
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh.thisDb(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE,
+            IOobject::NO_REGISTER
+        ),
+        mesh
+    );
+
+    // Expresions of volFields
+    {
+        volScalarField result
+        (
+            IOobject
+            (
+                "result",
+                runTime.timeName(),
+                mesh.thisDb(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                IOobject::NO_REGISTER
+            ),
+            mesh,
+            dimensionedScalar(p.dimensions(), 0)
+        );
+
+        {
+            Pout<< "No expression templates:" << endl;
+            const high_resolution_clock::time_point t1 =
+                high_resolution_clock::now();
+            result = p + p2;
+            const high_resolution_clock::time_point t2 =
+                high_resolution_clock::now();
+            const duration<double> time_span = t2 - t1;
+            Pout<< "Operation time:" << time_span.count() << endl;
+        }
+        {
+            Pout<< "With expression templates:" << endl;
+            const high_resolution_clock::time_point t1 =
+                high_resolution_clock::now();
+            result = p.expr() + p2.expr();
+            const high_resolution_clock::time_point t2 =
+                high_resolution_clock::now();
+            const duration<double> time_span = t2 - t1;
+            Pout<< "Operation time:" << time_span.count() << endl;
+        }
+
+//        const auto oldDimensions = p.dimensions();
+//        p.dimensions().reset(dimless);
+//        p2.dimensions().reset(dimless);
+//        result.dimensions().reset(dimless);
+//        {
+//
+//            Pout<< "Complex expression : No expression templates:" << endl;
+//            const high_resolution_clock::time_point t1 =
+//                high_resolution_clock::now();
+//            result = cos(p + 0.5*sqrt(p2-sin(p)));
+//            const high_resolution_clock::time_point t2 =
+//                high_resolution_clock::now();
+//            const duration<double> time_span = t2 - t1;
+//            Pout<< "Operation time:" << time_span.count() << endl;
+//        }
+//        {
+//            Pout<< "Complex expression : With expression templates:" << endl;
+//            const high_resolution_clock::time_point t1 =
+//                high_resolution_clock::now();
+//            const auto zeroDotFive
+//            (
+//                dimensionedScalar(dimless, 0.5).expr(p)
+//            );
+//            result = cos(p.expr() + zeroDotFive*sqrt(p2.expr()-sin(p.expr())));
+//            const high_resolution_clock::time_point t2 =
+//                high_resolution_clock::now();
+//            const duration<double> time_span = t2 - t1;
+//            Pout<< "Operation time:" << time_span.count() << endl;
+//        }
+//        p.dimensions().reset(oldDimensions);
+//        p2.dimensions().reset(oldDimensions);
+//        result.dimensions().reset(oldDimensions);
+
+        return 0;
+//        auto expression = someFunction(p).expr() + someFunction(p).expr();
+//        result = expression;
+//        DebugVar(result);
+    }
+    {
+        volScalarField result
+        (
+            IOobject
+            (
+                "result",
+                runTime.timeName(),
+                mesh.thisDb(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                IOobject::NO_REGISTER
+            ),
+            mesh,
+            dimensionedScalar(p.dimensions(), 0)
+        );
+        auto expression = someFunction(p).expr() + someFunction(p).expr();
+        result = expression;
+        DebugVar(result);
+    }
+
+    // Expresions of volFields
+    {
+        volScalarField result
+        (
+            "result",
+            mesh,
+            sqr(p.expr() + p.expr())
+        );
+        DebugVar(result);
+    }
+    {
+        // Fill p with some values
+        forAll(p, celli)
+        {
+            p[celli] = celli;
+        }
+        p.correctBoundaryConditions();
+
+        // Interpolate to surface field
+        surfaceScalarField result
+        (
+            "result",
+            mesh,
+            Expression::interpolate    //<volExpr, surfaceExpr>
+            (
+                p.expr(),
+                mesh.surfaceInterpolation::weights().expr(),
+                mesh
+            )
+        );
+        DebugVar(result);
+    }
+    {
+        // For testing as a replacement of laplacian weights
+        const volScalarField gamma
+        (
+            IOobject
+            (
+                "gamma",
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                IOobject::NO_REGISTER
+            ),
+            mesh,
+            dimensionedScalar(dimless, 1.0)
+        );
+
+        fvMatrix<scalar> fvm
+        (
+            p,
+            gamma.dimensions()*mesh.magSf().dimensions()*p.dimensions()
+        );
+
+        const linear<scalar> interpGammaScheme(mesh);
+        const fv::correctedSnGrad<scalar> snGradScheme(mesh);
+
+        fusedGaussFvmLaplacian
+        (
+            fvm,
+            interpGammaScheme,
+            snGradScheme,
+            gamma,
+            p
+        );
+
+        DebugVar(fvm.source());
+    }
+
+
+
+    // Expressions of fvMatrix
+    {
+        tmp<fvMatrix<scalar>> tm0(fvm::laplacian(p));
+        const fvMatrix<scalar>& m0 = tm0();
+        DebugVar(m0.dimensions());
+
+        tmp<fvMatrix<scalar>> tm1(fvm::laplacian(p));
+        const fvMatrix<scalar>& m1 = tm1();
+        DebugVar(m1.dimensions());
+
+        fvMatrix<scalar> m2(p, m0.expr() + m1.expr());
+        DebugVar(m2.dimensions());
+    }
+
+    return 0;
+}
+
+
+// ************************************************************************* //

src/Allwmake

@@ -55,8 +55,6 @@ wmake $targetType fileFormats
 wmake $targetType surfMesh
 wmake $targetType meshTools
 
-#------------------------------------------------------------------------------
-
 wmake $targetType finiteArea
 wmake $targetType finiteVolume
 wmake $targetType fused/finiteVolume

src/OpenFOAM/containers/Lists/List/List.H

@@ -45,6 +45,7 @@ SourceFiles
 
 #include "autoPtr.H"
 #include "UList.H"
+#include "SLListFwd.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -65,7 +66,6 @@ typedef List<bool> boolList;        //!< A List of bools
 typedef List<char> charList;        //!< A List of chars
 typedef List<label> labelList;      //!< A List of labels
 
-
 /*---------------------------------------------------------------------------*\
                            Class List Declaration
 \*---------------------------------------------------------------------------*/
@@ -395,6 +395,37 @@ public:
 
         //- Alias for resize()
         void setSize(label n, const T& val) { this->resize(n, val); }
+
+
+    // Expression templates
+
+        //- Construct from value expression
+        template<typename E>
+        explicit List
+        (
+            const Expression::ListExpression
+            <
+                E,
+                typename E::value_type
+            >& expr
+        )
+        {
+            expr.evaluate(*this);
+        }
+
+        //- Assign values from expression
+        template<typename E>
+        void operator=
+        (
+            const Expression::ListExpression
+            <
+                E,
+                typename E::value_type
+            >& expr
+        )
+        {
+            expr.evaluate(*this);
+        }
 };
 
 

src/OpenFOAM/containers/Lists/List/UList.H

@@ -80,6 +80,11 @@ typedef UList<bool> boolUList;      //!< A UList of bools
 typedef UList<char> charUList;      //!< A UList of chars
 typedef UList<label> labelUList;    //!< A UList of labels
 
+namespace Expression
+{
+template<class E, typename ValType> class ListExpression;
+template<class T> class ListConstRefWrap;
+}
 
 /*---------------------------------------------------------------------------*\
                            Class UList Declaration
@@ -697,6 +702,39 @@ public:
         {
             return this->contains(val, pos);
         }
+
+
+    // Expression templates
+
+        //- Have unique tag
+        using is_List = void;
+
+        //- Wrap value as expression
+        auto expr() const
+        {
+            return Expression::ListConstRefWrap<T>
+            (
+                reinterpret_cast<const List<T>&>(*this)
+            );
+        }
+
+        //- Assign values from expression
+        template<typename E>
+        void operator=
+        (
+            const Expression::ListExpression
+            <
+                E,
+                typename E::value_type
+            >& expr
+        )
+        {
+            // bypass expr.evaluate to avoid resize ...
+            for (label i = 0; i < expr.size(); ++i)
+            {
+                operator[](i) = expr[i];
+            }
+        }
 };
 
 

src/OpenFOAM/dimensionedTypes/dimensionedType/dimensionedType.H

@@ -60,6 +60,12 @@ template<class Type> class dimensioned;
 template<class Type>
 Istream& operator>>(Istream& is, dimensioned<Type>& dt);
 
+namespace Expression
+{
+template<class T> class UniformListWrap;
+template<class GeoField> class UniformGeometricFieldWrap;
+}
+
 /*---------------------------------------------------------------------------*\
                            Class dimensioned Declaration
 \*---------------------------------------------------------------------------*/
@@ -449,6 +455,28 @@ public:
         {
             return getOrAddToDict(name, dict, deflt);
         }
+
+
+    // Expression templates
+
+        //- Wrap value as constant-value List expression
+        auto expr(const label size) const
+        {
+            return Expression::UniformListWrap<Type>(size, value_);
+        }
+
+        //- Wrap value as constant-value GeometricField expression. Supplied
+        //  fields is used for sizing only and needs to be valid at time of
+        //  evaluation.
+        template<class GeoField>
+        auto expr(const GeoField& fld) const
+        {
+            return Expression::UniformGeometricFieldWrap<GeoField>
+            (
+                fld,
+                *this
+            );
+        }
 };
 
 

src/OpenFOAM/expressionTemplates/GenericExpression.H

@@ -0,0 +1,252 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2025 M. Janssens
+-------------------------------------------------------------------------------
+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/>.
+
+Namespace
+    Foam::GenericExpression
+
+Description
+    Expression templates.
+
+SourceFiles
+    GenericExpression.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_GenericExpression_H
+#define Foam_GenericExpression_H
+
+#include <cassert>
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace Expression
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class GenericExpression Declarations
+\*---------------------------------------------------------------------------*/
+
+template<typename E>
+class GenericExpression
+{
+public:
+    static constexpr bool is_leaf = false;
+
+    auto evaluate() const
+    {
+        // Delegation to the actual expression type. This avoids dynamic
+        // polymorphism (a.k.a. virtual functions in C++)
+        return static_cast<E const&>(*this).evaluate();
+    }
+};
+
+
+// Macros
+// ~~~~~~
+
+#define EXPRESSION_FUNCTION1(BaseClass, Func, BaseFunc, OpFunc)                \
+template<typename E1>                                                          \
+class OpFunc                                                                   \
+:                                                                              \
+    public BaseClass<OpFunc<E1>>                                               \
+{                                                                              \
+    typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_;      \
+                                                                               \
+public:                                                                        \
+    static constexpr bool is_leaf = false;                                     \
+                                                                               \
+    OpFunc(const E1& u)                                                        \
+    :                                                                          \
+        u_(u)                                                                  \
+    {}                                                                         \
+                                                                               \
+    auto evaluate() const                                                      \
+    {                                                                          \
+        return BaseFunc(u_.evaluate());                                        \
+    }                                                                          \
+};                                                                             \
+template<typename E1>                                                          \
+OpFunc<E1> Func                                                                \
+(                                                                              \
+    const BaseClass<E1>& u                                                     \
+)                                                                              \
+{                                                                              \
+    return OpFunc<E1>(static_cast<const E1&>(u));                              \
+}
+
+
+#define EXPRESSION_FUNCTION2(BaseClass, Func, BaseFunc, OpFunc)                \
+template<typename E1, typename E2>                                             \
+class OpFunc                                                                   \
+:                                                                              \
+    public BaseClass<OpFunc<E1, E2>>                                           \
+{                                                                              \
+    typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_;      \
+    typename std::conditional<E2::is_leaf, const E2&, const E2>::type v_;      \
+                                                                               \
+public:                                                                        \
+    static constexpr bool is_leaf = false;                                     \
+                                                                               \
+    OpFunc(const E1& u, const E2& v)                                           \
+    :                                                                          \
+     u_(u), v_(v)                                                              \
+    {}                                                                         \
+    auto evaluate() const                                                      \
+    {                                                                          \
+        return BaseFunc(u_.evaluate(), v_.evaluate());                         \
+    }                                                                          \
+};                                                                             \
+template<typename E1, typename E2>                                             \
+OpFunc<E1, E2> Func                                                            \
+(                                                                              \
+    const BaseClass<E1>& u,                                                    \
+    const BaseClass<E2>& v                                                     \
+)                                                                              \
+{                                                                              \
+    return OpFunc<E1, E2>                                                      \
+    (                                                                          \
+        static_cast<const E1&>(u),                                             \
+        static_cast<const E2&>(v)                                              \
+    );                                                                         \
+}
+
+
+#define EXPRESSION_OPERATOR(BaseClass, Op, BaseOp, OpFunc)                     \
+template<typename E1, typename E2>                                             \
+class OpFunc                                                                   \
+:                                                                              \
+    public BaseClass                                                           \
+    <                                                                          \
+        OpFunc<E1, E2>                                                         \
+    >                                                                          \
+{                                                                              \
+    /* cref if leaf, copy otherwise */                                         \
+    typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_;      \
+    typename std::conditional<E2::is_leaf, const E2&, const E2>::type v_;      \
+                                                                               \
+public:                                                                        \
+    static constexpr bool is_leaf = false;                                     \
+                                                                               \
+    OpFunc(const E1& u, const E2& v)                                           \
+    :                                                                          \
+     u_(u), v_(v)                                                              \
+    {}                                                                         \
+    auto evaluate() const                                                      \
+    {                                                                          \
+        return u_.evaluate() BaseOp v_.evaluate();                             \
+    }                                                                          \
+};                                                                             \
+template<typename E1, typename E2>                                             \
+OpFunc<E1, E2>                                                                 \
+operator Op                                                                    \
+(                                                                              \
+    BaseClass<E1> const& u,                                                    \
+    BaseClass<E2> const& v                                                     \
+)                                                                              \
+{                                                                              \
+    return OpFunc<E1, E2>                                                      \
+    (                                                                          \
+        static_cast<const E1&>(u),                                             \
+        static_cast<const E2&>(v)                                              \
+    );                                                                         \
+}
+
+
+// Do '-' separately until we work out macro expansion...
+template<typename E1>
+class g_negate
+:
+    public GenericExpression
+    <
+        g_negate<E1>
+    >
+{
+    typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_;
+
+public:
+    static constexpr bool is_leaf = false;
+
+    g_negate(const E1& u)
+    :
+        u_(u)
+    {}
+
+    auto evaluate() const
+    {
+        return -u_.evaluate();
+    }
+};
+template<typename E1>
+g_negate<E1> operator-
+(
+    const GenericExpression<E1>& u
+)
+{
+    return g_negate<E1>(static_cast<const E1&>(u));
+}
+
+
+EXPRESSION_FUNCTION1(GenericExpression, sqr, Foam::sqr, g_sqr)
+EXPRESSION_FUNCTION1(GenericExpression, sqrt, Foam::sqrt, g_sqrt)
+EXPRESSION_FUNCTION1(GenericExpression, magSqr, Foam::magSqr, g_magSqr)
+EXPRESSION_FUNCTION1(GenericExpression, symm, Foam::symm, g_symm)
+EXPRESSION_FUNCTION1(GenericExpression, pow2, Foam::pow2, g_pow2)
+EXPRESSION_FUNCTION1(GenericExpression, pow3, Foam::pow3, g_pow3)
+EXPRESSION_FUNCTION1(GenericExpression, pow4, Foam::pow4, g_pow4)
+EXPRESSION_FUNCTION1(GenericExpression, operator~, Foam::operator~, g_compl)
+
+//TBD. Parse problem
+//EXPRESSION_FUNCTION1(GenericExpression, operator!, Foam::operator!, g_not)
+
+#undef EXPRESSION_FUNCTION1
+
+EXPRESSION_FUNCTION2(GenericExpression, operator+, Foam::operator+, g_add);
+EXPRESSION_FUNCTION2(GenericExpression, operator-, Foam::operator-, g_subtract);
+EXPRESSION_FUNCTION2(GenericExpression, operator*, Foam::operator*, g_multiply);
+EXPRESSION_FUNCTION2(GenericExpression, operator/, Foam::operator/, g_divide);
+
+#undef EXPRESSION_FUNCTION2
+
+EXPRESSION_OPERATOR(GenericExpression, ||, ||, g_or)
+EXPRESSION_OPERATOR(GenericExpression, &&, &&, g_and)
+EXPRESSION_OPERATOR(GenericExpression, &, &, g_bitand)
+EXPRESSION_OPERATOR(GenericExpression, |, |, g_bitor)
+EXPRESSION_OPERATOR(GenericExpression, ^, ^, g_xor)
+
+#undef EXPRESSION_OPERATOR
+
+} // End namespace Expression
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/expressionTemplates/boolExpression.H

@@ -0,0 +1,145 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2025 M. Janssens
+-------------------------------------------------------------------------------
+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/>.
+
+Namespace
+    Foam::boolExpression
+
+Description
+    Expression templates for bool
+
+SourceFiles
+    boolExpression.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_boolExpression_H
+#define Foam_boolExpression_H
+
+#include "GenericExpression.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace Expression
+{
+
+/*---------------------------------------------------------------------------*\
+                       Namespace Expression Declarations
+\*---------------------------------------------------------------------------*/
+
+// Wrap of non-const bool
+// ~~~~~~~~~~~~~~~~~~~~~~
+
+class boolWrap
+:
+    public GenericExpression<boolWrap>
+{
+    bool elems_;
+
+public:
+    static constexpr bool is_leaf = false;  //true;
+
+    //- Construct from components
+    boolWrap(const bool elems)
+    :
+        elems_(elems)
+    {}
+
+    // Construct from GenericExpression, forcing its evaluation.
+    template<typename E>
+    boolWrap
+    (
+        const GenericExpression<E>& expr
+    )
+    :
+        elems_(expr.evaluate())
+    {}
+
+    auto evaluate() const
+    {
+        return elems_;
+    }
+
+    template<typename E>
+    auto evaluate
+    (
+        const GenericExpression<E>& expr
+    )
+    {
+        elems_ = expr.evaluate();
+        return elems_;
+    }
+
+    //- Assignment
+    template<typename E>
+    void operator=
+    (
+        const GenericExpression<E>& expr
+    )
+    {
+        elems_ = expr.evaluate();
+    }
+};
+
+
+// Wrap of const bool
+// ~~~~~~~~~~~~~~~~~~
+
+class boolConstWrap
+:
+    public GenericExpression<boolConstWrap>
+{
+    const bool elems_;
+
+public:
+
+    // ! Store as copy (since holds reference)
+    static constexpr bool is_leaf = false;
+
+    // construct from components
+    boolConstWrap(const bool elems)
+    :
+        elems_(elems)
+    {}
+
+    auto evaluate() const
+    {
+        return elems_;
+    }
+};
+
+
+} // End namespace Expression
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/expressionTemplates/dimensionSetExpression.H

@@ -0,0 +1,148 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2025 M. Janssens
+-------------------------------------------------------------------------------
+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/>.
+
+Namespace
+    Foam::dimensionSetExpression
+
+Description
+    Expression templates for dimensionSet
+
+SourceFiles
+    dimensionSetExpression.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_dimensionSetExpression_H
+#define Foam_dimensionSetExpression_H
+
+#include "GenericExpression.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace Expression
+{
+
+/*---------------------------------------------------------------------------*\
+                       Namespace Expression Declarations
+\*---------------------------------------------------------------------------*/
+
+// Wrap of non-const reference to dimensionSet
+// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+class dimensionSetRefWrap
+:
+    public GenericExpression<dimensionSetRefWrap>
+{
+    dimensionSet& elems_;
+
+public:
+    static constexpr bool is_leaf = false;  //true;
+
+    //- Construct from components
+    dimensionSetRefWrap(dimensionSet& elems)
+    :
+        elems_(elems)
+    {}
+
+    // Construct from GenericExpression, forcing its evaluation.
+    template<typename E>
+    dimensionSetRefWrap
+    (
+        dimensionSet& elems,
+        const GenericExpression<E>& expr
+    )
+    :
+        elems_(elems)
+    {
+        elems_ = expr.evaluate();
+    }
+
+    auto evaluate() const
+    {
+        return elems_;
+    }
+
+    template<typename E>
+    auto& evaluate
+    (
+        const GenericExpression<E>& expr
+    )
+    {
+        elems_ = expr.evaluate();
+        return elems_;
+    }
+
+    //- Assignment
+    template<typename E>
+    void operator=
+    (
+        const GenericExpression<E>& expr
+    )
+    {
+        elems_ = expr.evaluate();
+    }
+};
+
+
+// Wrap of const reference to dimensionSet
+// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+class dimensionSetConstRefWrap
+:
+    public GenericExpression<dimensionSetConstRefWrap>
+{
+    const dimensionSet& elems_;
+
+public:
+
+    // ! Store as copy (since holds reference)
+    static constexpr bool is_leaf = false;
+
+    // construct from components
+    dimensionSetConstRefWrap(const dimensionSet& elems)
+    :
+        elems_(elems)
+    {}
+
+    auto evaluate() const
+    {
+        return elems_;
+    }
+};
+
+
+} // End namespace Expression
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/fields/Fields/Field/Field.H

@@ -582,6 +582,37 @@ public:
 
         friend Ostream& operator<< <Type>
         (Ostream&, const tmp<Field<Type>>&);
+
+
+    // Expression templates
+
+        //- Construct from value expression
+        template<typename E>
+        explicit Field
+        (
+            const Expression::ListExpression
+            <
+                E,
+                typename E::value_type
+            >& expr
+        )
+        {
+            expr.evaluate(*this);
+        }
+
+        //- Assign values from expression
+        template<typename E>
+        void operator=
+        (
+            const Expression::ListExpression
+            <
+                E,
+                typename E::value_type
+            >& expr
+        )
+        {
+            expr.evaluate(*this);
+        }
 };
 
 

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C

@@ -1675,6 +1675,124 @@ Foam::Ostream& Foam::operator<<
 }
 
 
+// * * * * * * * * * * * * * Expression Templates  * * * * * * * * * * * * * //
+
+template<class Type, template<class> class PatchField, class GeoMesh>
+template<typename E>
+Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
+(
+    const IOobject& io,
+    const Mesh& mesh,
+    const Expression::GeometricFieldExpression
+    <
+        E,
+        typename E::IntExpr,
+        typename E::UncoupledPatchExpr,
+        typename E::CoupledPatchExpr,
+        typename E::value_type
+    >& expr
+)
+:
+    Internal(io, mesh, expr.dimensions(), false),
+    timeIndex_(this->time().timeIndex()),
+    boundaryField_(mesh.boundary(), *this, PatchField<Type>::calculatedType())
+{
+    DebugInFunction
+        << "Creating from expression " << nl << this->info() << endl;
+
+    bool hasRead = readIfPresent();
+    if (!hasRead)
+    {
+        expr.evaluate(*this);
+    }
+}
+
+
+template<class Type, template<class> class PatchField, class GeoMesh>
+template<typename E>
+Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
+(
+    const word& name,
+    const Mesh& mesh,
+    const Expression::GeometricFieldExpression
+    <
+        E,
+        typename E::IntExpr,
+        typename E::UncoupledPatchExpr,
+        typename E::CoupledPatchExpr,
+        typename E::value_type
+    >& expr
+)
+:
+    Internal
+    (
+        IOobject
+        (
+            name,
+            mesh.time().timeName(),
+            mesh.thisDb()
+        ),
+        mesh,
+        expr.dimensions(),
+        false
+    ),
+    timeIndex_(this->time().timeIndex()),
+    boundaryField_(mesh.boundary(), *this, PatchField<Type>::calculatedType())
+{
+    DebugInFunction
+        << "Creating from expression " << nl << this->info() << endl;
+
+    expr.evaluate(*this);
+}
+
+
+template<class Type, template<class> class PatchField, class GeoMesh>
+Foam::Expression::GeometricFieldConstRefWrap
+<
+    Foam::GeometricField<Type, PatchField, GeoMesh>
+>
+Foam::GeometricField<Type, PatchField, GeoMesh>::expr() const
+{
+    return Expression::GeometricFieldConstRefWrap<this_type>(*this);
+}
+
+
+template<class Type, template<class> class PatchField, class GeoMesh>
+template<typename E>
+void Foam::GeometricField<Type, PatchField, GeoMesh>::operator=
+(
+    const Expression::GeometricFieldExpression
+    <
+        E,
+        typename E::IntExpr,
+        typename E::UncoupledPatchExpr,
+        typename E::CoupledPatchExpr,
+        typename E::value_type
+    >& expr
+)
+{
+    Expression::GeometricFieldRefWrap<this_type>(*this, expr);
+}
+
+
+template<class Type, template<class> class PatchField, class GeoMesh>
+template<typename E>
+void Foam::GeometricField<Type, PatchField, GeoMesh>::operator==
+(
+    const Expression::GeometricFieldExpression
+    <
+        E,
+        typename E::IntExpr,
+        typename E::UncoupledPatchExpr,
+        typename E::CoupledPatchExpr,
+        typename E::value_type
+    >& expr
+)
+{
+    expr.evaluate(*this, true);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #undef checkField

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H

@@ -42,6 +42,7 @@ SourceFiles
 #define Foam_GeometricField_H
 
 #include "GeometricBoundaryField.H"
+#include "GeometricFieldExpression.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -1142,6 +1143,79 @@ public:
         {
             return this->clamp_range(lo.value(), hi.value());
         }
+
+
+    // Expression templates
+
+        //- Have unique tag
+        using is_GeometricField = void;
+
+        //- Construct from IOobject and value expression. Supports
+        //  read_if_present
+        template<typename E>
+        GeometricField
+        (
+            const IOobject& io,
+            const Mesh& mesh,
+            const Expression::GeometricFieldExpression
+            <
+                E,
+                typename E::IntExpr,
+                typename E::UncoupledPatchExpr,
+                typename E::CoupledPatchExpr,
+                typename E::value_type
+            >& expr
+        );
+
+        //- Construct from name, mesh and value expression.
+        template<typename E>
+        GeometricField
+        (
+            const word& name,
+            const Mesh& mesh,
+            const Expression::GeometricFieldExpression
+            <
+                E,
+                typename E::IntExpr,
+                typename E::UncoupledPatchExpr,
+                typename E::CoupledPatchExpr,
+                typename E::value_type
+            >& expr
+        );
+
+        //- Wrap value as expression
+        typename Expression::GeometricFieldConstRefWrap
+        <
+            this_type
+        > expr() const;
+
+        //- Assign values from expression
+        template<typename E>
+        void operator=
+        (
+            const Expression::GeometricFieldExpression
+            <
+                E,
+                typename E::IntExpr,
+                typename E::UncoupledPatchExpr,
+                typename E::CoupledPatchExpr,
+                typename E::value_type
+            >& expr
+        );
+
+        //- Assign values from expression. Override asssignable.
+        template<typename E>
+        void operator==
+        (
+            const Expression::GeometricFieldExpression
+            <
+                E,
+                typename E::IntExpr,
+                typename E::UncoupledPatchExpr,
+                typename E::CoupledPatchExpr,
+                typename E::value_type
+            >& expr
+        );
 };
 
 

src/OpenFOAM/memory/tmp/tmp.H

@@ -59,8 +59,17 @@ See also
 namespace Foam
 {
 
+namespace Expression
+{
+template<class GeoField> class GeometricFieldConstTmpWrap;
+template<class T> class ListConstTmpWrap;
+}
+
 // Forward Declarations
 template<class T> class refPtr;
+template<class T> class tmp;
+template<class Type, template<class> class PatchField, class GeoMesh>
+class GeometricField;
 
 /*---------------------------------------------------------------------------*\
                              Class tmp Declaration
@@ -352,6 +361,33 @@ public:
         //  \deprecated(2020-07) - use bool operator
         FOAM_DEPRECATED_FOR(2020-07, "bool operator")
         bool empty() const noexcept { return !ptr_; }
+
+
+    // Expression templates
+
+        //- Wrap value as expression
+        auto expr() const
+        {
+            // Add 'using is_GeometricField = void' to GeometricField
+            // Add 'using is_List = void' to List
+
+            if constexpr
+            (
+                std::is_void_v<typename T::is_GeometricField>
+            )
+            {
+                return Expression::GeometricFieldConstTmpWrap<T>(*this);
+            }
+            else //if constexpr (std::is_void_v<T::is_List>)
+            {
+                return Expression::ListConstTmpWrap<T>(*this);
+            }
+            //else
+            //{
+            //    //static_assert("XXX",YYY);
+            //    return;
+            //}
+        }
 };
 
 

src/OpenFOAM/memory/tmp/tmpI.H

@@ -131,7 +131,7 @@ inline Foam::tmp<T>::tmp(const tmp<T>& rhs)
         if (ptr_)
         {
             ptr_->refCount::operator++();
-            this->checkUseCount();
+            //this->checkUseCount();
         }
         else
         {
@@ -162,7 +162,7 @@ inline Foam::tmp<T>::tmp(const tmp<T>& rhs, bool reuse)
             else
             {
                 ptr_->refCount::operator++();
-                this->checkUseCount();
+                //this->checkUseCount();
             }
         }
         else

src/finiteVolume/finiteVolume/gradSchemes/LeastSquaresGrad/LeastSquaresGrad.C

@@ -1,4 +1,4 @@
-/*---------------------------------------------------------------------------*\
+  /*---------------------------------------------------------------------------*\
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
@@ -36,23 +36,18 @@ License
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 template<class Type, class Stencil>
-Foam::tmp
-<
-    Foam::GeometricField
-    <
-        typename Foam::outerProduct<Foam::vector, Type>::type,
-        Foam::fvPatchField,
-    Foam::volMesh
-    >
->
-Foam::fv::LeastSquaresGrad<Type, Stencil>::calcGrad
+void Foam::fv::LeastSquaresGrad<Type, Stencil>::calcGrad
 (
-    const GeometricField<Type, fvPatchField, volMesh>& vtf,
-    const word& name
+    GeometricField
+    <
+        typename outerProduct<vector, Type>::type,
+        fvPatchField,
+        volMesh
+    >& lsGrad,
+    const GeometricField<Type, fvPatchField, volMesh>& vtf
 ) const
 {
     typedef typename outerProduct<vector, Type>::type GradType;
-    typedef GeometricField<GradType, fvPatchField, volMesh> GradFieldType;
 
     const fvMesh& mesh = vtf.mesh();
 
@@ -62,24 +57,7 @@ Foam::fv::LeastSquaresGrad<Type, Stencil>::calcGrad
         mesh
     );
 
-    tmp<GradFieldType> tlsGrad
-    (
-        new GradFieldType
-        (
-            IOobject
-            (
-                name,
-                vtf.instance(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh,
-            dimensioned<GradType>(vtf.dimensions()/dimLength, Zero),
-            fvPatchFieldBase::extrapolatedCalculatedType()
-        )
-    );
-    GradFieldType& lsGrad = tlsGrad.ref();
+    lsGrad = dimensioned<GradType>(vtf.dimensions()/dimLength, Zero);
     Field<GradType>& lsGradIf = lsGrad;
 
     const extendedCentredCellToCellStencil& stencil = lsv.stencil();
@@ -131,6 +109,50 @@ Foam::fv::LeastSquaresGrad<Type, Stencil>::calcGrad
     // Correct the boundary conditions
     lsGrad.correctBoundaryConditions();
     gaussGrad<Type>::correctBoundaryConditions(vtf, lsGrad);
+}
+
+
+template<class Type, class Stencil>
+Foam::tmp
+<
+    Foam::GeometricField
+    <
+        typename Foam::outerProduct<Foam::vector, Type>::type,
+        Foam::fvPatchField,
+    Foam::volMesh
+    >
+>
+Foam::fv::LeastSquaresGrad<Type, Stencil>::calcGrad
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vtf,
+    const word& name
+) const
+{
+    typedef typename outerProduct<vector, Type>::type GradType;
+    typedef GeometricField<GradType, fvPatchField, volMesh> GradFieldType;
+
+    const fvMesh& mesh = vtf.mesh();
+
+    tmp<GradFieldType> tlsGrad
+    (
+        new GradFieldType
+        (
+            IOobject
+            (
+                name,
+                vtf.instance(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            vtf.dimensions()/dimLength,
+            fvPatchFieldBase::extrapolatedCalculatedType()
+        )
+    );
+    GradFieldType& lsGrad = tlsGrad.ref();
+
+    calcGrad(lsGrad, vtf);
 
     return tlsGrad;
 }

src/finiteVolume/finiteVolume/gradSchemes/LeastSquaresGrad/LeastSquaresGrad.H

@@ -132,6 +132,14 @@ public:
             const GeometricField<Type, fvPatchField, volMesh>& vsf,
             const word& name
         ) const;
+
+        //- Calculate the grad of the given field into supplied field
+        virtual void calcGrad
+        (
+            GeometricField
+            <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>& res,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const;
 };
 
 

src/finiteVolume/finiteVolume/gradSchemes/fourthGrad/fourthGrad.C

@@ -38,19 +38,21 @@ License
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 template<class Type>
-Foam::tmp
-<
-    Foam::GeometricField
-    <
-        typename Foam::outerProduct<Foam::vector, Type>::type,
-        Foam::fvPatchField,
-        Foam::volMesh
-    >
->
-Foam::fv::fourthGrad<Type>::calcGrad
+void Foam::fv::fourthGrad<Type>::calcGrad
 (
+    GeometricField
+    <
+        typename outerProduct<vector, Type>::type,
+        fvPatchField,
+        volMesh
+    >& fGrad,
     const GeometricField<Type, fvPatchField, volMesh>& vsf,
-    const word& name
+    const GeometricField
+    <
+        typename outerProduct<vector, Type>::type,
+        fvPatchField,
+        volMesh
+    >& secondfGrad
 ) const
 {
     // The fourth-order gradient is calculated in two passes.  First,
@@ -59,37 +61,11 @@ Foam::fv::fourthGrad<Type>::calcGrad
     // gradient to complete the accuracy.
 
     typedef typename outerProduct<vector, Type>::type GradType;
-    typedef GeometricField<GradType, fvPatchField, volMesh> GradFieldType;
 
     const fvMesh& mesh = vsf.mesh();
 
-    // Assemble the second-order least-square gradient
-    // Calculate the second-order least-square gradient
-    tmp<GradFieldType> tsecondfGrad
-      = leastSquaresGrad<Type>(mesh).grad
-        (
-            vsf,
-            "leastSquaresGrad(" + vsf.name() + ")"
-        );
-    const GradFieldType& secondfGrad =
-        tsecondfGrad();
-
-    tmp<GradFieldType> tfGrad
-    (
-        new GradFieldType
-        (
-            IOobject
-            (
-                name,
-                vsf.instance(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            secondfGrad
-        )
-    );
-    GradFieldType& fGrad = tfGrad.ref();
+    // Initialise to gradient
+    fGrad = secondfGrad;
 
     const vectorField& C = mesh.C();
 
@@ -158,9 +134,102 @@ Foam::fv::fourthGrad<Type>::calcGrad
 
     fGrad.correctBoundaryConditions();
     gaussGrad<Type>::correctBoundaryConditions(vsf, fGrad);
+}
+
+
+template<class Type>
+Foam::tmp
+<
+    Foam::GeometricField
+    <
+        typename Foam::outerProduct<Foam::vector, Type>::type,
+        Foam::fvPatchField,
+        Foam::volMesh
+    >
+>
+Foam::fv::fourthGrad<Type>::calcGrad
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vsf,
+    const word& name
+) const
+{
+    // The fourth-order gradient is calculated in two passes.  First,
+    // the standard least-square gradient is assembled.  Then, the
+    // fourth-order correction is added to the second-order accurate
+    // gradient to complete the accuracy.
+
+    typedef typename outerProduct<vector, Type>::type GradType;
+    typedef GeometricField<GradType, fvPatchField, volMesh> GradFieldType;
+
+    const fvMesh& mesh = vsf.mesh();
+
+    // Assemble the second-order least-square gradient
+    // Calculate the second-order least-square gradient
+    tmp<GradFieldType> tsecondfGrad
+      = leastSquaresGrad<Type>(mesh).grad
+        (
+            vsf,
+            "leastSquaresGrad(" + vsf.name() + ")"
+        );
+    const GradFieldType& secondfGrad = tsecondfGrad();
+
+    tmp<GradFieldType> tfGrad
+    (
+        new GradFieldType
+        (
+            IOobject
+            (
+                name,
+                vsf.instance(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            vsf.dimensions()/dimLength,
+            fvPatchFieldBase::extrapolatedCalculatedType()
+        )
+    );
+
+    calcGrad(tfGrad.ref(), vsf, secondfGrad);
 
     return tfGrad;
 }
 
 
+template<class Type>
+void Foam::fv::fourthGrad<Type>::calcGrad
+(
+    GeometricField
+    <
+        typename outerProduct<vector, Type>::type,
+        fvPatchField,
+        volMesh
+    >& fGrad,
+    const GeometricField<Type, fvPatchField, volMesh>& vsf
+) const
+{
+    // The fourth-order gradient is calculated in two passes.  First,
+    // the standard least-square gradient is assembled.  Then, the
+    // fourth-order correction is added to the second-order accurate
+    // gradient to complete the accuracy.
+
+    typedef typename outerProduct<vector, Type>::type GradType;
+    typedef GeometricField<GradType, fvPatchField, volMesh> GradFieldType;
+
+    const fvMesh& mesh = vsf.mesh();
+
+    // Assemble the second-order least-square gradient
+    // Calculate the second-order least-square gradient
+    tmp<GradFieldType> tsecondfGrad
+      = leastSquaresGrad<Type>(mesh).grad
+        (
+            vsf,
+            "leastSquaresGrad(" + vsf.name() + ")"
+        );
+
+    calcGrad(fGrad, vsf, tsecondfGrad());
+}
+
+
 // ************************************************************************* //

src/finiteVolume/finiteVolume/gradSchemes/fourthGrad/fourthGrad.H

@@ -99,6 +99,32 @@ public:
             const GeometricField<Type, fvPatchField, volMesh>& vsf,
             const word& name
         ) const;
+
+        //- Calculate the grad of the given field into supplied field
+        virtual void calcGrad
+        (
+            GeometricField
+            <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>& res,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const;
+
+        //- Helper : calculate the grad of the given field into supplied field
+        void calcGrad
+        (
+            GeometricField
+            <
+                typename outerProduct<vector, Type>::type,
+                fvPatchField,
+                volMesh
+            >& fGrad,
+            const GeometricField<Type, fvPatchField, volMesh>& vsf,
+            const GeometricField
+            <
+                typename outerProduct<vector, Type>::type,
+                fvPatchField,
+                volMesh
+            >& secondfGrad
+        ) const;
 };
 
 

src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C

@@ -140,6 +140,67 @@ Foam::fv::gaussGrad<Type>::calcGrad
 }
 
 
+template<class Type>
+void Foam::fv::gaussGrad<Type>::calcGrad
+(
+    GeometricField
+    <
+        typename outerProduct<vector, Type>::type,
+        fvPatchField,
+        volMesh
+    >& gGrad,
+    const GeometricField<Type, fvPatchField, volMesh>& vsf
+) const
+{
+    typedef typename outerProduct<vector, Type>::type GradType;
+
+    DebugPout<< "gaussGrad<Type>::calcGrad on " << vsf.name()
+        << " into " << gGrad.name() << endl;
+
+    const fvMesh& mesh = vsf.mesh();
+    const labelUList& owner = mesh.owner();
+    const labelUList& neighbour = mesh.neighbour();
+    const vectorField& Sf = mesh.Sf();
+
+    const auto tssf(tinterpScheme_().interpolate(vsf));
+    const auto& ssf = tssf();
+
+    gGrad = dimensioned<GradType>(vsf.dimensions()/dimLength, Zero);
+
+    Field<GradType>& igGrad = gGrad;
+    const Field<Type>& issf = ssf;
+
+    forAll(owner, facei)
+    {
+        const GradType Sfssf = Sf[facei]*issf[facei];
+
+        igGrad[owner[facei]] += Sfssf;
+        igGrad[neighbour[facei]] -= Sfssf;
+    }
+
+    forAll(mesh.boundary(), patchi)
+    {
+        const labelUList& pFaceCells =
+            mesh.boundary()[patchi].faceCells();
+
+        const vectorField& pSf = mesh.Sf().boundaryField()[patchi];
+
+        const fvsPatchField<Type>& pssf = ssf.boundaryField()[patchi];
+
+        forAll(mesh.boundary()[patchi], facei)
+        {
+            igGrad[pFaceCells[facei]] += pSf[facei]*pssf[facei];
+        }
+    }
+
+    igGrad /= mesh.V();
+
+    gGrad.correctBoundaryConditions();
+
+    correctBoundaryConditions(vsf, gGrad);
+}
+
+
 template<class Type>
 void Foam::fv::gaussGrad<Type>::correctBoundaryConditions
 (

src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.H

@@ -145,6 +145,14 @@ public:
             const word& name
         ) const;
 
+        //- Calculate the grad of the given field into supplied field
+        virtual void calcGrad
+        (
+            GeometricField
+            <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>& res,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const;
+
         //- Correct the boundary values of the gradient using the patchField
         //- snGrad functions
         static void correctBoundaryConditions

src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.C

@@ -96,7 +96,7 @@ Foam::fv::gradScheme<Type>::grad
     GradFieldType* pgGrad =
         mesh().objectRegistry::template getObjectPtr<GradFieldType>(name);
 
-    if (!this->mesh().cache(name) || this->mesh().changing())
+    if (!this->mesh().cache(name) || this->mesh().topoChanging())
     {
         // Delete any old occurrences to avoid double registration
         if (pgGrad && pgGrad->ownedByRegistry())
@@ -119,17 +119,14 @@ Foam::fv::gradScheme<Type>::grad
     }
     else
     {
-        if (pgGrad->upToDate(vsf))
+        if (pgGrad->upToDate(vsf) && this->mesh().upToDatePoints(*pgGrad))
         {
             solution::cachePrintMessage("Reusing", name, vsf);
         }
         else
         {
             solution::cachePrintMessage("Updating", name, vsf);
-            delete pgGrad;
-
-            pgGrad = calcGrad(vsf, name).ptr();
-            regIOobject::store(pgGrad);
+            calcGrad(*pgGrad, vsf);
         }
     }
 
@@ -180,4 +177,23 @@ Foam::fv::gradScheme<Type>::grad
 }
 
 
+template<class Type>
+void Foam::fv::gradScheme<Type>::calcGrad
+(
+    GeometricField
+    <
+        typename outerProduct<vector, Type>::type,
+        fvPatchField,
+        volMesh
+    >& result,
+    const GeometricField<Type, fvPatchField, volMesh>& vsf
+) const
+{
+    DebugPout<< "gradScheme<Type>::calcGrad on " << vsf.name()
+        << " into " << result.name() << endl;
+
+    result = calcGrad(vsf, result.name());
+}
+
+
 // ************************************************************************* //

src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.H

@@ -142,6 +142,14 @@ public:
             const word& name
         ) const = 0;
 
+        //- Calculate the grad of the given field into supplied field
+        virtual void calcGrad
+        (
+            GeometricField
+            <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>& res,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const;
+
         //- Calculate and return the grad of the given field
         //- which may have been cached
         tmp

src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrad.C

@@ -62,40 +62,25 @@ void Foam::fv::cellLimitedGrad<Type, Limiter>::limitGradient
 
 
 template<class Type, class Limiter>
-Foam::tmp
-<
-    Foam::GeometricField
-    <
-        typename Foam::outerProduct<Foam::vector, Type>::type,
-        Foam::fvPatchField,
-        Foam::volMesh
-    >
->
-Foam::fv::cellLimitedGrad<Type, Limiter>::calcGrad
+void Foam::fv::cellLimitedGrad<Type, Limiter>::calcGrad
 (
-    const GeometricField<Type, fvPatchField, volMesh>& vsf,
-    const word& name
-) const
-{
-    const fvMesh& mesh = vsf.mesh();
-
-    tmp
-    <
-        GeometricField
-        <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
-    > tGrad = basicGradScheme_().calcGrad(vsf, name);
-
-    if (k_ < SMALL)
-    {
-        return tGrad;
-    }
-
     GeometricField
     <
         typename outerProduct<vector, Type>::type,
         fvPatchField,
         volMesh
-    >& g = tGrad.ref();
+    >& g,
+    const GeometricField<Type, fvPatchField, volMesh>& vsf
+) const
+{
+    const fvMesh& mesh = vsf.mesh();
+
+    basicGradScheme_().calcGrad(g, vsf);
+
+    if (k_ < SMALL)
+    {
+        return;
+    }
 
     const labelUList& owner = mesh.owner();
     const labelUList& neighbour = mesh.neighbour();
@@ -227,6 +212,49 @@ Foam::fv::cellLimitedGrad<Type, Limiter>::calcGrad
     limitGradient(limiter, g);
     g.correctBoundaryConditions();
     gaussGrad<Type>::correctBoundaryConditions(vsf, g);
+}
+
+
+template<class Type, class Limiter>
+Foam::tmp
+<
+    Foam::GeometricField
+    <
+        typename Foam::outerProduct<Foam::vector, Type>::type,
+        Foam::fvPatchField,
+        Foam::volMesh
+    >
+>
+Foam::fv::cellLimitedGrad<Type, Limiter>::calcGrad
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vsf,
+    const word& name
+) const
+{
+    typedef typename outerProduct<vector, Type>::type GradType;
+    typedef GeometricField<GradType, fvPatchField, volMesh> GradFieldType;
+
+    const fvMesh& mesh = vsf.mesh();
+
+    tmp<GradFieldType> tGrad
+    (
+        new GradFieldType
+        (
+            IOobject
+            (
+                name,
+                vsf.instance(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            vsf.dimensions()/dimLength,
+            fvPatchFieldBase::extrapolatedCalculatedType()
+        )
+    );
+    
+    calcGrad(tGrad.ref(), vsf);
 
     return tGrad;
 }

src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/cellLimitedGrad.H

@@ -150,6 +150,14 @@ public:
             const GeometricField<Type, fvPatchField, volMesh>& vsf,
             const word& name
         ) const;
+
+        //- Calculate the grad of the given field into supplied field
+        virtual void calcGrad
+        (
+            GeometricField
+            <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>& res,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const;
 };
 
 

src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C

@@ -2984,6 +2984,77 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const fvMatrix<Type>& fvm)
 }
 
 
+// * * * * * * * * * * * * * Expression Templates  * * * * * * * * * * * * * //
+
+template<class Type>
+template<typename E>
+Foam::fvMatrix<Type>::fvMatrix
+(
+    const GeometricField<Type, fvPatchField, volMesh>& psi,
+    const Expression::fvMatrixExpression
+    <
+        E,
+        typename E::DiagExpr,
+        typename E::UpperExpr,
+        typename E::LowerExpr,
+        typename E::FaceFluxExpr,
+        typename E::SourceExpr
+   >& expr
+)
+:
+    lduMatrix(psi.mesh()),
+    psi_(psi),
+    useImplicit_(false),
+    lduAssemblyName_(),
+    nMatrix_(0),
+    dimensions_(expr.dimensions()),
+    source_(psi.size(), Zero),
+    internalCoeffs_(psi.mesh().boundary().size()),
+    boundaryCoeffs_(psi.mesh().boundary().size())
+{
+    DebugInFunction
+        << "Constructing fvMatrix<Type> from expression for field "
+        << psi_.name() << endl;
+
+    checkImplicit();
+
+    // Fill in diag,upper,lower etc.
+    expr.evaluate(*this);
+
+    auto& psiRef = this->psi(0);
+    const label currentStatePsi = psiRef.eventNo();
+    psiRef.boundaryFieldRef().updateCoeffs();
+    psiRef.eventNo() = currentStatePsi;
+}
+
+
+template<class Type>
+Foam::Expression::fvMatrixConstRefWrap<Foam::fvMatrix<Type>>
+Foam::fvMatrix<Type>::expr() const
+{
+    return Expression::fvMatrixConstRefWrap<Foam::fvMatrix<Type>>(*this);
+}
+
+
+template<class Type>
+template<typename E>
+void Foam::fvMatrix<Type>::operator=
+(
+    const Expression::fvMatrixExpression
+    <
+        E,
+        typename E::DiagExpr,
+        typename E::UpperExpr,
+        typename E::LowerExpr,
+        typename E::FaceFluxExpr,
+        typename E::SourceExpr
+   >& expr
+)
+{
+    expr.evaluate(*this);
+}
+
+
 // * * * * * * * * * * * * * * * * Solvers * * * * * * * * * * * * * * * * * //
 
 #include "fvMatrixSolve.C"

src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H

@@ -53,6 +53,8 @@ SourceFiles
 #include "lduPrimitiveMeshAssembly.H"
 #include "lduMesh.H"
 
+#include "fvMatrixExpression.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -741,6 +743,47 @@ public:
             Ostream&,
             const fvMatrix<Type>&
         );
+
+
+    // Expression templates
+
+        //- Construct given a field to solve for and expressions for
+        //- all the components
+        template<typename E>
+        fvMatrix
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& psi,
+            const Expression::fvMatrixExpression
+            <
+                E,
+                typename E::DiagExpr,
+                typename E::UpperExpr,
+                typename E::LowerExpr,
+                typename E::FaceFluxExpr,
+                typename E::SourceExpr
+           >& expr
+        );
+
+        //- Wrap value as expression
+        typename Expression::fvMatrixConstRefWrap
+        <
+            fvMatrix<Type>
+        > expr() const;
+
+        //- Assign values from expression
+        template<typename E>
+        void operator=
+        (
+            const Expression::fvMatrixExpression
+            <
+                E,
+                typename E::DiagExpr,
+                typename E::UpperExpr,
+                typename E::LowerExpr,
+                typename E::FaceFluxExpr,
+                typename E::SourceExpr
+           >& expr
+       );
 };
 
 

src/functionObjects/field/Make/files

@@ -108,10 +108,6 @@ turbulenceFields/turbulenceFields.C
 yPlus/yPlus.C
 wallShearStress/wallShearStress.C
 wallHeatFlux/wallHeatFlux.C
-wallHeatFlux/wallHeatFluxModels/wallHeatFluxModel/wallHeatFluxModel.C
-wallHeatFlux/wallHeatFluxModels/wallHeatFluxModel/wallHeatFluxModelNew.C
-wallHeatFlux/wallHeatFluxModels/wall/wall.C
-wallHeatFlux/wallHeatFluxModels/gauge/gauge.C
 
 writeCellCentres/writeCellCentres.C
 writeCellVolumes/writeCellVolumes.C

src/functionObjects/field/wallHeatFlux/wallHeatFlux.C

@@ -5,7 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2025 OpenCFD Ltd.
+    Copyright (C) 2016-2017 OpenFOAM Foundation
+    Copyright (C) 2016-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -26,8 +27,14 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "wallHeatFlux.H"
-#include "wallHeatFluxModel.H"
+#include "turbulentFluidThermoModel.H"
+#include "solidThermo.H"
+#include "surfaceInterpolate.H"
+#include "fvcSnGrad.H"
+#include "wallPolyPatch.H"
+#include "turbulentFluidThermoModel.H"
 #include "addToRunTimeSelectionTable.H"
+#include "multiphaseInterSystem.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -40,6 +47,54 @@ namespace functionObjects
 }
 }
 
+
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+void Foam::functionObjects::wallHeatFlux::writeFileHeader(Ostream& os) const
+{
+    writeHeader(os, "Wall heat-flux");
+    writeCommented(os, "Time");
+    writeTabbed(os, "patch");
+    writeTabbed(os, "min");
+    writeTabbed(os, "max");
+    writeTabbed(os, "integral");
+    os  << endl;
+}
+
+
+void Foam::functionObjects::wallHeatFlux::calcHeatFlux
+(
+    const volScalarField& alpha,
+    const volScalarField& he,
+    volScalarField& wallHeatFlux
+)
+{
+    volScalarField::Boundary& wallHeatFluxBf = wallHeatFlux.boundaryFieldRef();
+
+    const volScalarField::Boundary& heBf = he.boundaryField();
+
+    const volScalarField::Boundary& alphaBf = alpha.boundaryField();
+
+    for (const label patchi : patchIDs_)
+    {
+        wallHeatFluxBf[patchi] = alphaBf[patchi]*heBf[patchi].snGrad();
+    }
+
+
+    const auto* qrPtr = cfindObject<volScalarField>(qrName_);
+
+    if (qrPtr)
+    {
+        const volScalarField::Boundary& radHeatFluxBf = qrPtr->boundaryField();
+
+        for (const label patchi : patchIDs_)
+        {
+            wallHeatFluxBf[patchi] -= radHeatFluxBf[patchi];
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::functionObjects::wallHeatFlux::wallHeatFlux
@@ -50,19 +105,32 @@ Foam::functionObjects::wallHeatFlux::wallHeatFlux
 )
 :
     fvMeshFunctionObject(name, runTime, dict),
-    qModelPtr_
-    (
-        wallHeatFluxModel::New
-        (
-            dict,
-            mesh_,
-            name,
-            scopedName(typeName),
-            *this
-        )
-    )
+    writeFile(obr_, name, typeName, dict),
+    qrName_("qr")
 {
     read(dict);
+
+    writeFileHeader(file());
+
+    volScalarField* wallHeatFluxPtr
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                scopedName(typeName),
+                mesh_.time().timeName(),
+                mesh_.thisDb(),
+                IOobjectOption::NO_READ,
+                IOobjectOption::NO_WRITE,
+                IOobjectOption::REGISTER
+            ),
+            mesh_,
+            dimensionedScalar(dimMass/pow3(dimTime), Zero)
+        )
+    );
+
+    regIOobject::store(wallHeatFluxPtr);
 }
 
 
@@ -70,26 +138,178 @@ Foam::functionObjects::wallHeatFlux::wallHeatFlux
 
 bool Foam::functionObjects::wallHeatFlux::read(const dictionary& dict)
 {
-    Log << type() << ' ' << name() << " read:" << endl;
-    if (!fvMeshFunctionObject::read(dict) || !qModelPtr_->read(dict))
+    const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
+
+    fvMeshFunctionObject::read(dict);
+    writeFile::read(dict);
+
+    dict.readIfPresent("qr", qrName_);
+
+    wordRes patchNames;
+    labelHashSet patchSet;
+    if (dict.readIfPresent("patches", patchNames) && !patchNames.empty())
     {
-        return false;
+        patchSet = pbm.patchSet(patchNames);
     }
+
+    labelHashSet allWalls(pbm.findPatchIDs<wallPolyPatch>());
+
+    Info<< type() << ' ' << name() << ':' << nl;
+
+    if (patchSet.empty())
+    {
+        patchIDs_ = allWalls.sortedToc();
+
+        Info<< "    processing all (" << patchIDs_.size()
+            << ") wall patches" << nl << endl;
+    }
+    else
+    {
+        allWalls &= patchSet;
+        patchSet -= allWalls;
+        patchIDs_ = allWalls.sortedToc();
+
+        if (!patchSet.empty())
+        {
+            WarningInFunction
+                << "Requested wall heat-flux on ("
+                << patchSet.size() << ") non-wall patches:" << nl;
+
+            for (const label patchi : patchSet.sortedToc())
+            {
+                Info<< "        " << pbm[patchi].name() << nl;
+            }
+            Info<< nl;
+        }
+
+        Info<< "    processing (" << patchIDs_.size()
+            << ") wall patches:" << nl;
+
+        for (const label patchi : patchIDs_)
+        {
+            Info<< "        " << pbm[patchi].name() << nl;
+        }
+        Info<< endl;
+    }
+
     return true;
 }
 
 
 bool Foam::functionObjects::wallHeatFlux::execute()
 {
-    Log << type() << ' ' << name() << " execute:" << endl;
-    return qModelPtr_->execute();
+    auto& wallHeatFlux = lookupObjectRef<volScalarField>(scopedName(typeName));
+
+    if
+    (
+        foundObject<compressible::turbulenceModel>
+        (
+            turbulenceModel::propertiesName
+        )
+    )
+    {
+        const compressible::turbulenceModel& turbModel =
+            lookupObject<compressible::turbulenceModel>
+            (
+                turbulenceModel::propertiesName
+            );
+
+        calcHeatFlux
+        (
+            turbModel.alphaEff()(),
+            turbModel.transport().he(),
+            wallHeatFlux
+        );
+    }
+    else if (foundObject<fluidThermo>(fluidThermo::dictName))
+    {
+        const fluidThermo& thermo =
+            lookupObject<fluidThermo>(fluidThermo::dictName);
+
+        calcHeatFlux
+        (
+            thermo.alpha(),
+            thermo.he(),
+            wallHeatFlux
+        );
+    }
+    else if (foundObject<solidThermo>(solidThermo::dictName))
+    {
+        const solidThermo& thermo =
+            lookupObject<solidThermo>(solidThermo::dictName);
+
+        calcHeatFlux(thermo.alpha(), thermo.he(), wallHeatFlux);
+    }
+    else if
+    (
+        foundObject<multiphaseInterSystem>
+            (multiphaseInterSystem::phasePropertiesName)
+    )
+    {
+        const auto& thermo = lookupObject<multiphaseInterSystem>
+        (
+            multiphaseInterSystem::phasePropertiesName
+        );
+
+        calcHeatFlux(thermo.kappaEff()(), thermo.T(), wallHeatFlux);
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Unable to find compressible turbulence model in the "
+            << "database" << exit(FatalError);
+    }
+
+
+    const fvPatchList& patches = mesh_.boundary();
+
+    const surfaceScalarField::Boundary& magSf = mesh_.magSf().boundaryField();
+
+    for (const label patchi : patchIDs_)
+    {
+        const fvPatch& pp = patches[patchi];
+
+        const scalarField& hfp = wallHeatFlux.boundaryField()[patchi];
+
+        const MinMax<scalar> limits = gMinMax(hfp);
+        const scalar integralHfp = gWeightedSum(magSf[patchi], hfp);
+
+        if (Pstream::master())
+        {
+            writeCurrentTime(file());
+
+            file()
+                << token::TAB << pp.name()
+                << token::TAB << limits.min()
+                << token::TAB << limits.max()
+                << token::TAB << integralHfp
+                << endl;
+        }
+
+        Log << "    min/max/integ(" << pp.name() << ") = "
+            << limits.min() << ", " << limits.max()
+            << ", " << integralHfp << endl;
+
+        this->setResult("min(" + pp.name() + ")", limits.min());
+        this->setResult("max(" + pp.name() + ")", limits.max());
+        this->setResult("int(" + pp.name() + ")", integralHfp);
+    }
+
+    return true;
 }
 
 
 bool Foam::functionObjects::wallHeatFlux::write()
 {
-    Log << type() << ' ' << name() << " write:" << endl;
-    return qModelPtr_->write();
+    const auto& wallHeatFlux =
+        lookupObject<volScalarField>(scopedName(typeName));
+
+    Log << type() << ' ' << name() << " write:" << nl
+        << "    writing field " << wallHeatFlux.name() << endl;
+
+    wallHeatFlux.write();
+
+    return true;
 }
 
 

src/functionObjects/field/wallHeatFlux/wallHeatFlux.H

@@ -5,7 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2025 OpenCFD Ltd.
+    Copyright (C) 2016-2017 OpenFOAM Foundation
+    Copyright (C) 2016-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,47 +28,31 @@ Class
     Foam::functionObjects::wallHeatFlux
 
 Group
+    grpForcesFunctionObjects
 
 Description
-    This function object calculates the wall heat flux based on a custom model
-    defined in the dictionary. The model can be specified to compute the heat
-    flux according to user-defined criteria or equations.
+    Computes the wall-heat flux at selected wall patches.
+
+    \table
+      Operand        | Type              | Location
+      input          | -                 | -
+      output file    | dat               | postProcessing/<FO>/<time>/field
+      output field   | volScalarField (only boundaryField) <!--
+                     --> | postProcessing/<FO>/<time>/outField
+    \endtable
 
 Usage
     Minimal example by using \c system/controlDict.functions:
     \verbatim
-    wallHeatFlux1
+    wallHeatFluxFO
     {
         // Mandatory entries
         type        wallHeatFlux;
         libs        (fieldFunctionObjects);
 
         // Optional entries
-        model       <word>;
-
-        // Conditional entries
-
-            // Option-1: when 'model' is 'wall' - previous 'wallHeatFlux' state
-
-            // Optional entries
-            patches     (<patch1> ... <patchN>); // (wall1 "(wall2|wall3)");
-            qr          <word>;
-
-
-            // Option-2: when 'model' is 'gauge'
-            // Mandatory entries
-            patch       <word>;
-            Tgauge      <scalar>;
-
-            // Optional entries
-            absorptivity  <scalar>;
-            emissivity    <scalar>;
-            T             <word>;
-            qin           <word>;
-            alphat        <word>;
-            convective    <bool>;
-            radiative     <bool>;
-            writeFields   <bool>;
+        patches     (<patch1> ... <patchN>); // (wall1 "(wall2|wall3)");
+        qr          <word>;;
 
         // Inherited entries
         ...
@@ -79,45 +64,31 @@ Usage
       Property     | Description                           | Type | Reqd | Deflt
       type         | Type name: wallHeatFlux               | word |  yes  | -
       libs         | Library name: fieldFunctionObjects    | word |  yes  | -
-      model      | Model name                              | word |  no  | wall
       patches    | Names of operand patches   | wordList | no | all wall patches
       qr         | Name of radiative heat flux field       | word | no | qr
-      patch      | Name of the patch to probe              | word | yes  | -
-      Tgauge     | Gauge temperature (K)                  | scalar | yes | -
-      absorptivity | Absorptivity of the gauge              | scalar | no | 1
-      emissivity   | Emissivity of the gauge                | scalar | no | 1
-      T           | Name of the temperature field          | word | no  | T
-      qin | Name of the incident radiative heat flux field | word | no | qin
-      alphat | Name of turbulent thermal diffusivity field | word | no | alphat
-      convective  | Write convective heat flux into file   | bool | no | true
-      radiative   | Write radiative heat flux into file    | bool | no | true
-      writeFields | Write the fields 'qConv' and 'qRad'    | bool | no | true
     \endtable
 
     The inherited entries are elaborated in:
-     - \link fvMeshFunctionObject.H \endlink
-     - \link writeFile.H \endlink
-     - \link pointProberBase.H \endlink
-     - \link patchProber.H \endlink
+      - \link functionObject.H \endlink
+      - \link writeFile.H \endlink
 
 SourceFiles
     wallHeatFlux.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef functionObjects_wallHeatFlux_H
-#define functionObjects_wallHeatFlux_H
+#ifndef Foam_functionObjects_wallHeatFlux_H
+#define Foam_functionObjects_wallHeatFlux_H
 
 #include "fvMeshFunctionObject.H"
+#include "writeFile.H"
+#include "volFieldsFwd.H"
+#include "HashSet.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
-// Forward Declarations
-class wallHeatFluxModel;
-
 namespace functionObjects
 {
 
@@ -127,12 +98,32 @@ namespace functionObjects
 
 class wallHeatFlux
 :
-    public fvMeshFunctionObject
+    public fvMeshFunctionObject,
+    public writeFile
 {
-    // Private Data
+protected:
 
-        //- Heat-flux model
-        autoPtr<wallHeatFluxModel> qModelPtr_;
+    // Protected Data
+
+        //- Wall patches to process (optionally filtered by name)
+        labelList patchIDs_;
+
+        //- Name of radiative heat flux name
+        word qrName_;
+
+
+    // Protected Member Functions
+
+        //- File header information
+        virtual void writeFileHeader(Ostream& os) const;
+
+        //- Calculate the heat-flux
+        void calcHeatFlux
+        (
+            const volScalarField& alpha,
+            const volScalarField& he,
+            volScalarField& wallHeatFlux
+        );
 
 
 public:
@@ -148,9 +139,15 @@ public:
         (
             const word& name,
             const Time& runTime,
-            const dictionary& dict
+            const dictionary&
         );
 
+        //- No copy construct
+        wallHeatFlux(const wallHeatFlux&) = delete;
+
+        //- No copy assignment
+        void operator=(const wallHeatFlux&) = delete;
+
 
     //- Destructor
     virtual ~wallHeatFlux() = default;
@@ -158,13 +155,13 @@ public:
 
     // Member Functions
 
-        //- Read the settings
+        //- Read the function-object dictionary
         virtual bool read(const dictionary& dict);
 
-        //- Calculate the heat-flux data
+        //- Execute the function-object operations
         virtual bool execute();
 
-        //- Write the heat-flux data
+        //- Write the function-object results
         virtual bool write();
 };
 

src/functionObjects/field/wallHeatFlux/wallHeatFluxModels/gauge/gauge.C

@@ -1,548 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2025 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 "gauge.H"
-#include "physicoChemicalConstants.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace wallHeatFluxModels
-{
-    defineTypeNameAndDebug(gauge, 0);
-    addToRunTimeSelectionTable
-    (
-        wallHeatFluxModel,
-        gauge,
-        dictionary
-    );
-}
-}
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-Foam::scalar Foam::wallHeatFluxModels::gauge::calcRadiantExitance() const
-{
-    // Radiant exitance (Wikipedia EN): https://w.wiki/CdZz
-    return e_*constant::physicoChemical::sigma.value()*Foam::pow4(Tgauge_);
-}
-
-
-Foam::tmp<Foam::labelField>
-Foam::wallHeatFluxModels::gauge::identifyProbeCells() const
-{
-    // Fetch the neighbour cells of the operand patch
-    const polyPatch& patch = mesh().boundaryMesh()[patchID_];
-    const labelUList& faceCells = patch.faceCells();
-    const label patchFaceStart = patch.start();
-
-    // Fetch the patch faces that correspond to the specified probes
-    const labelList& probeFaces = patchProber::faces();
-
-    // Initialise the size of the operand cells
-    auto tprobeCells = tmp<labelField>::New(probeFaces.size(), -1);
-    labelField& probeCells = tprobeCells.ref();
-
-    // Store the indices of the cells that contain the probed patch faces
-    forAll(probeFaces, facei)
-    {
-        if (activeFaces_[facei])
-        {
-            const label patchFaceGlobali = probeFaces[facei];
-            const label patchFaceLocali = patchFaceGlobali - patchFaceStart;
-
-            probeCells[facei] = faceCells[patchFaceLocali];
-        }
-    }
-
-    return tprobeCells;
-}
-
-
-bool Foam::wallHeatFluxModels::gauge::writeFileHeader(Ostream& os)
-{
-    writeHeader(os, typeName);
-    os  << "# Patch," << mesh().boundaryMesh()[patchID_].name() << nl
-        << "# Tgauge," << Tgauge_ << nl
-        << "# Absorptivity," << a_ << nl
-        << "# Emissivity," << e_ << nl
-        << "# ProbeID,Original Location,Patch Location" << nl;
-
-    const pointField& oldPoints = patchProber::oldPoints();
-    const pointField& probeLocs = patchProber::probeLocations();
-
-    for (label i = 0; i < szProbes_; ++i)
-    {
-        const vector& oldPoint = oldPoints[i];
-        const vector& probeLoc = probeLocs[i];
-
-        os  << '#' << ' ' << i
-            << ',' << oldPoint.x() << ',' << oldPoint.y() << ',' << oldPoint.z()
-            << ',' << probeLoc.x() << ',' << probeLoc.y() << ',' << probeLoc.z()
-            << nl;
-    }
-
-    os  << "# Time";
-    for (int i = 0; i < szProbes_; ++i)
-    {
-        os  << ',' << i;
-    }
-
-    os  << endl;
-    return true;
-}
-
-
-bool Foam::wallHeatFluxModels::gauge::calcConvectiveHeatFlux()
-{
-    // Fetch the turbulent thermal diffusivity field
-    const auto* alphatPtr = mesh().cfindObject<volScalarField>(alphatName_);
-
-#ifdef FULLDEBUG
-    if (!alphatPtr)
-    {
-        FatalErrorInFunction
-            << "Turbulent thermal diffusivity field, alphat, is not available."
-            << nl << "alphat: " << alphatName_
-            << exit(FatalError);
-    }
-#endif
-
-    // Compute the term: C_p*(alpha + alpha_t)
-    tmp<volScalarField> tkappaEff = thermo_.kappaEff(*alphatPtr);
-    const volScalarField& kappaEff = tkappaEff.cref();
-
-    // Sample the term above at patch-probe locations
-    tmp<scalarField> tsampledKappaEff = patchProber::sample(kappaEff);
-    const scalarField& sampledKappaEff = tsampledKappaEff.cref();
-
-
-    // Compute the patch-normal gradient of temperature with respect to the
-    // specified gauge temperature
-    tmp<scalarField> tdTdn = calcdTdn();
-    const scalarField& dTdn = tdTdn.cref();
-
-#ifdef FULLDEBUG
-    if (sampledKappaEff.size() != dTdn.size())
-    {
-        FatalErrorInFunction
-            << "Size mismatch: kappaEff samples (" << sampledKappaEff.size()
-            << ") vs dTdn (" << dTdn.size() << ")." << nl
-            << "Probe selection/order must match." << nl
-            << exit(FatalError);
-    }
-#endif
-
-    // Compute q''_conv = k_eff * dT/dn
-    tmp<scalarField> tq = sampledKappaEff*dTdn;
-    const scalarField& q = tq.cref();
-
-    // Be picky - Ensure storage exists and matches size
-    if (qConvPtr_ && (qConvPtr_->size() == q.size()))
-    {
-        *qConvPtr_ = q;
-    }
-    else
-    {
-        IOobject io
-        (
-            "qConv",
-            mesh().time().timeName(),
-            mesh()
-        );
-
-        qConvPtr_.reset(new scalarIOField(io, q));
-    }
-
-
-    // Note that for inactive faces q=0 since dT/dn=0
-    Pstream::listReduce(*qConvPtr_, sumOp<scalar>());
-
-    return true;
-}
-
-
-Foam::tmp<Foam::scalarField>
-Foam::wallHeatFluxModels::gauge::calcdTdn() const
-{
-    // Allocate and initialise the storage for dT/dn
-    auto tdTdn = tmp<scalarField>::New(szProbes_, 0);
-    scalarField& dTdn = tdTdn.ref();
-
-    // Fetch the temperature field
-    const auto* Tptr = mesh().cfindObject<volScalarField>(Tname_);
-
-#ifdef FULLDEBUG
-    if (!Tptr)
-    {
-        FatalErrorInFunction
-            << "Temperature field, T, is not available."
-            << nl << "T: " << Tname_
-            << exit(FatalError);
-    }
-#endif
-
-    // Fetch the internal fields for temperature and distance-to-patch fields
-    const scalarField& Ti = Tptr->primitiveField();
-    const auto& pdc = mesh().deltaCoeffs().boundaryField()[patchID_];
-
-
-    // Compute dT/dn for active faces and their owner cells
-    forAll(cells_, probei)
-    {
-        if (activeFaces_[probei])
-        {
-            const label cellID = cells_[probei];
-            dTdn[probei] = (Ti[cellID] - Tgauge_)/max(pdc[probei], SMALL);
-        }
-    }
-
-    return tdTdn;
-}
-
-
-bool Foam::wallHeatFluxModels::gauge::calcRadiativeHeatFlux()
-{
-    // Fetch the incident radiative heat-flux field
-    const auto* qinPtr = mesh().getObjectPtr<volScalarField>(qinName_);
-
-#ifdef FULLDEBUG
-    if (!qinPtr)
-    {
-        FatalErrorInFunction
-            << "Incident radiative heat-flux field, qin, is not available."
-            << nl << "qin: " << qinName_
-            << exit(FatalError);
-    }
-#endif
-
-    const volScalarField& qin = *qinPtr;
-
-    // Sample the incident radiative heat-flux field at patch-probe locations
-    tmp<scalarField> tq  = patchProber::sample(qin);
-    const scalarField& q = tq.cref();
-
-    // Compute q''_rad = a * q_in - M_e
-    tmp<scalarField> tqrad = a_*q - Me_;
-    scalarField& qrad = tqrad.ref();
-
-    // Overwrite q''_rad=0 for inactive probed faces
-    forAll(qrad, probei)
-    {
-        if (!activeFaces_[probei])
-        {
-            qrad[probei] = 0;
-        }
-    }
-
-
-    // Be picky - Ensure storage exists and matches size
-    if (qRadPtr_ && (qRadPtr_->size() == qrad.size()))
-    {
-        *qRadPtr_ = qrad;
-    }
-    else
-    {
-        IOobject io
-        (
-            "qRad",
-            mesh().time().timeName(),
-            mesh()
-        );
-
-        qRadPtr_.reset(new scalarIOField(io, qrad));
-    }
-
-    // Note that for inactive faces qrad=0
-    Pstream::listReduce(*qRadPtr_, sumOp<scalar>());
-
-    return true;
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::wallHeatFluxModels::gauge::gauge
-(
-    const dictionary& dict,
-    const fvMesh& mesh,
-    const word& name,
-    const word objName,
-    functionObjects::stateFunctionObject& state
-)
-:
-    wallHeatFluxModel(dict, mesh, name, objName, state),
-    patchProber(mesh, dict),
-    thermo_(mesh.lookupObject<fluidThermo>(fluidThermo::dictName))
-{}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-bool Foam::wallHeatFluxModels::gauge::read(const dictionary& dict)
-{
-    if (!wallHeatFluxModel::read(dict) || !patchProber::read(dict))
-    {
-        return false;
-    }
-
-
-    // Do not proceed if the operand heat-flux field is not registered
-    qinName_ = dict.getOrDefault<word>("qin", "qin");
-
-    const auto* qinPtr = mesh().cfindObject<volScalarField>(qinName_);
-    if (!qinPtr)
-    {
-        FatalErrorInFunction
-            << "Incident radiative heat-flux field, qin, is not available."
-            << nl << "qin: " << qinName_
-            << exit(FatalError);
-    }
-
-
-    // Fetch the patch-prober data
-    const labelList& patchIDs = patchProber::patchIDs();
-    if (patchIDs.size() == 1)
-    {
-        patchID_ = patchIDs[0];
-    }
-    else
-    {
-        FatalIOErrorInFunction(dict)
-            << "This function object can operate on only a single patch." << nl
-            << "The number of input patches: " << patchIDs.size()
-            << exit(FatalIOError);
-    }
-
-
-    // Fetch the number of probes, and verify the consistent sizing
-    szProbes_ = patchProber::probeLocations().size();
-
-    const labelList& probeFaces = patchProber::faces();
-
-    if (szProbes_ != probeFaces.size())
-    {
-        FatalErrorInFunction
-            << "Size mismatch: szProbes (" << szProbes_
-            << ") vs probeFaces (" << probeFaces.size() << ")." << nl
-            << "Number of patch probes and patch faces must match." << nl
-            << exit(FatalError);
-    }
-
-
-    // Create and register the output heat-flux fields
-    // Maybe number of probes is changed; thus, no check for the pointer
-    qConvPtr_.reset
-    (
-        new scalarIOField
-        (
-            IOobject
-            (
-                "qConv",
-                mesh().time().timeName(),
-                mesh()
-            ),
-            szProbes_
-        )
-    );
-
-    qRadPtr_.reset
-    (
-        new scalarIOField
-        (
-            IOobject
-            (
-                "qRad",
-                mesh().time().timeName(),
-                mesh()
-            ),
-            szProbes_
-        )
-    );
-
-
-    // Read and store the common properties of the gauges
-    Tgauge_ = dict.getScalar("Tgauge");
-    if (Tgauge_ < 0)
-    {
-        FatalIOErrorInFunction(dict)
-            << "Gauge temperature cannot be negative." << nl
-            << "Tgauge: " << Tgauge_
-            << exit(FatalIOError);
-    }
-
-    a_ = dict.getOrDefault<scalar>("absorptivity", 1);
-    e_ = dict.getOrDefault<scalar>("emissivity", 1);
-    if (a_ < 0 || a_ > 1 || e_ < 0 || e_ > 1)
-    {
-        FatalIOErrorInFunction(dict)
-            << "The range of absorptivity and emissivity can only be [0,1]."
-            << "Absorptivity: " << a_ << nl
-            << "Emissivity: " << e_
-            << exit(FatalIOError);
-    }
-
-
-    // Read names of various fields; also early check if the fields exist
-    Tname_ = dict.getOrDefault<word>("T", "T");
-    const auto* Tptr = mesh().cfindObject<volScalarField>(Tname_);
-    if (!Tptr)
-    {
-        FatalErrorInFunction
-            << "Temperature field, T, is not available."
-            << nl << "T: " << Tname_
-            << exit(FatalError);
-    }
-
-    alphatName_ = dict.getOrDefault<word>("alphat", "alphat");
-    const auto* alphatPtr = mesh().cfindObject<volScalarField>(alphatName_);
-    if (!alphatPtr)
-    {
-        FatalErrorInFunction
-            << "Turbulent thermal diffusivity field, alphat, is not available."
-            << nl << "alphat: " << alphatName_
-            << exit(FatalError);
-    }
-
-
-    // Compute the radiant-exitance term of the radiative heat flux
-    // The term is calculated once per simulation unless Tgauge is changed
-    Me_ = calcRadiantExitance();
-
-
-    // Identify the faces that are actively sampled - assuming no duplication
-    activeFaces_.resize_fill(probeFaces.size(), false);
-
-    forAll(probeFaces, facei)
-    {
-        if (probeFaces[facei] != -1)
-        {
-            activeFaces_[facei] = true;
-        }
-    }
-
-
-    // Identify the indices of the cells that contain the probed patch faces
-    cells_ = identifyProbeCells();
-
-
-    // Prepare the output files
-    if (writeFile::canResetFile())
-    {
-        writeFile::resetFile(objName());
-    }
-
-    if (writeFile::canWriteHeader())
-    {
-        writeFileHeader(file());
-        writtenHeader_ = true;
-    }
-
-    bool convective = dict.getOrDefault<bool>("convective", false);
-    if (UPstream::master() && convective && !convectiveFilePtr_)
-    {
-        convectiveFilePtr_ = newFileAtStartTime("convective");
-        writeFileHeader(convectiveFilePtr_());
-    }
-
-    bool radiative = dict.getOrDefault<bool>("radiative", false);
-    if (UPstream::master() && radiative && !radiativeFilePtr_)
-    {
-        radiativeFilePtr_ = newFileAtStartTime("radiative");
-        writeFileHeader(radiativeFilePtr_());
-    }
-
-
-    writeFields_ = dict.getOrDefault<bool>("writeFields", false);
-
-
-    return true;
-}
-
-
-bool Foam::wallHeatFluxModels::gauge::execute()
-{
-    if (!(calcConvectiveHeatFlux() && calcRadiativeHeatFlux()))
-    {
-        return false;
-    }
-    return true;
-}
-
-
-bool Foam::wallHeatFluxModels::gauge::write()
-{
-    if (UPstream::master())
-    {
-        Ostream& os = file();
-
-        os  << mesh().time().timeOutputValue();
-        for (label pi = 0; pi < szProbes_; ++pi)
-        {
-            os  << ',' << (*qConvPtr_)[pi] + (*qRadPtr_)[pi];
-        }
-        os  << endl;
-    }
-
-    if (UPstream::master() && convectiveFilePtr_)
-    {
-        OFstream& os = convectiveFilePtr_.ref();
-
-        os  << mesh().time().timeOutputValue();
-        for (label pi = 0; pi < szProbes_; ++pi)
-        {
-            os  << ',' << (*qConvPtr_)[pi];
-        }
-        os  << endl;
-    }
-
-    if (UPstream::master() && radiativeFilePtr_)
-    {
-        OFstream& os = radiativeFilePtr_.ref();
-
-        os  << mesh().time().timeOutputValue();
-        for (label pi = 0; pi < szProbes_; ++pi)
-        {
-            os  << ',' << (*qRadPtr_)[pi];
-        }
-        os  << endl;
-    }
-
-    if (writeFields_)
-    {
-        if (qConvPtr_) qConvPtr_->write();
-        if (qRadPtr_) qRadPtr_->write();
-    }
-
-    return true;
-}
-
-
-// ************************************************************************* //

src/functionObjects/field/wallHeatFlux/wallHeatFluxModels/gauge/gauge.H

@@ -1,247 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2025 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::wallHeatFluxModels::gauge
-
-Description
-    This model computes the wall-heat flux at a selected patch using
-    imaginary heat-flux gauges.
-
-    Heat-flux gauges are measurement tools that are placed on walls. The
-    temperature of the gauges is usually different from that of the surrounding
-    wall - an order of magnitude lower. The model allows the user to
-    specify the gauge temperature and other gauge properties such as
-    absorptivity or emissivity without creating separate patches.
-    The model calculates the net convective and radiative heat fluxes at the
-    gauge locations, combines them, and writes the results to the output files.
-
-Usage
-    Minimal example by using \c system/controlDict.functions:
-    \verbatim
-    wallHeatFlux1
-    {
-        // Mandatory entries
-        type        wallHeatFlux;
-        libs        (fieldFunctionObjects);
-
-        model       gauge;
-
-        Tgauge      <scalar>;
-        patch       <word>;
-
-        // Optional entries
-        absorptivity  <scalar>;
-        emissivity    <scalar>;
-        T             <word>;
-        qin           <word>;
-        alphat        <word>;
-        convective    <bool>;
-        radiative     <bool>;
-        writeFields   <bool>;
-
-        // Inherited entries
-        ...
-    }
-    \endverbatim
-
-    where the entries mean:
-    \table
-      Property     | Description                           | Type | Reqd | Deflt
-      type         | Type name: wallHeatFlux               | word |  yes  | -
-      libs         | Library name: fieldFunctionObjects    | word |  yes  | -
-      model      | Model name: gauge                     | word |  no  | wall
-      Tgauge     | Gauge temperature (K)                  | scalar | yes | -
-      patch      | Name of the patch to probe              | word | yes  | -
-      absorptivity | Absorptivity of the gauge              | scalar | no | 1
-      emissivity   | Emissivity of the gauge                | scalar | no | 1
-      T           | Name of the temperature field          | word | no  | T
-      qin | Name of the incident radiative heat flux field | word | no | qin
-      alphat | Name of turbulent thermal diffusivity field | word | no | alphat
-      convective  | Calculate convective heat flux         | bool | no | false
-      radiative   | Calculate radiative heat flux          | bool | no | false
-      writeFields | Write the fields to file                | bool | no | false
-    \endtable
-
-    The inherited entries are elaborated in:
-     - \link fvMeshFunctionObject.H \endlink
-     - \link writeFile.H \endlink
-     - \link pointProberBase.H \endlink
-     - \link patchProber.H \endlink
-
-SourceFiles
-    gauge.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef wallHeatFluxModels_gauge_H
-#define wallHeatFluxModels_gauge_H
-
-#include "wallHeatFluxModel.H"
-#include "writeFile.H"
-#include "patchProber.H"
-#include "wallDist.H"
-#include "turbulentFluidThermoModel.H"
-#include "scalarIOField.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace wallHeatFluxModels
-{
-
-/*---------------------------------------------------------------------------*\
-                            Class gauge Declaration
-\*---------------------------------------------------------------------------*/
-
-class gauge
-:
-    public wallHeatFluxModel,
-    public patchProber
-{
-    // Private Data
-
-        //- Const reference to the thermo database
-        const fluidThermo& thermo_;
-
-        //- Net convective gauge heat flux at patch-probe locations
-        std::unique_ptr<scalarIOField> qConvPtr_;
-
-        //- Net radiative gauge heat flux at patch-probe locations
-        std::unique_ptr<scalarIOField> qRadPtr_;
-
-        //- Flag to output the net convective heat flux data
-        autoPtr<OFstream> convectiveFilePtr_;
-
-        //- Flag to output the net radiative heat flux data
-        autoPtr<OFstream> radiativeFilePtr_;
-
-        //- Indices of the patch cells that will be probed
-        labelField cells_;
-
-        //- Identifier of the faces that are actively sampled
-        boolList activeFaces_;
-
-        //- Name of the incident radiative heat-flux field
-        word qinName_;
-
-        //- Name of the temperature field
-        word Tname_;
-
-        //- Name of the turbulent thermal diffusivity field
-        word alphatName_;
-
-        //- Common and constant temperature of the gauges
-        scalar Tgauge_;
-
-        //- Radiant exitance
-        scalar Me_;
-
-        //- Absorptivity of the gauge surfaces
-        scalar a_;
-
-        //- Emissivity of the gauge surfaces
-        scalar e_;
-
-        //- Index of the operand patch
-        label patchID_;
-
-        //- Number of specified patch probes
-        label szProbes_;
-
-        //- Flag to output the net convective/radiative heat flux fields
-        bool writeFields_;
-
-
-    // Private Member Functions
-
-        //- Return the radiant-exitance term of the radiative heat flux
-        scalar calcRadiantExitance() const;
-
-        //- Return the indices of the cells that contain the probed patch faces
-        tmp<labelField> identifyProbeCells() const;
-
-        //- Write the file-header information
-        bool writeFileHeader(Ostream& os);
-
-        //- Calculate the net convective heat flux
-        bool calcConvectiveHeatFlux();
-
-        //- Return the dT/dn term of the convective heat flux
-        tmp<scalarField> calcdTdn() const;
-
-        //- Calculate the net radiative heat flux
-        bool calcRadiativeHeatFlux();
-
-
-public:
-
-    //- Runtime type information
-    TypeName("gauge");
-
-
-    // Constructors
-
-        //- Construct from components
-        gauge
-        (
-            const dictionary& dict,
-            const fvMesh& mesh,
-            const word& name,
-            const word objName,
-            functionObjects::stateFunctionObject& state
-        );
-
-
-    //- Destructor
-    virtual ~gauge() = default;
-
-
-    // Member Functions
-
-    // Evaluation
-
-        //- Read the settings
-        virtual bool read(const dictionary& dict);
-
-        //- Calculate the heat-flux data
-        virtual bool execute();
-
-        //- Write the heat-flux data
-        virtual bool write();
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace wallHeatFluxModels
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/functionObjects/field/wallHeatFlux/wallHeatFluxModels/wall/wall.C

@@ -1,333 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2025 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 "wall.H"
-#include "turbulentFluidThermoModel.H"
-#include "solidThermo.H"
-#include "surfaceInterpolate.H"
-#include "fvcSnGrad.H"
-#include "wallPolyPatch.H"
-#include "turbulentFluidThermoModel.H"
-#include "multiphaseInterSystem.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace wallHeatFluxModels
-{
-    defineTypeNameAndDebug(wall, 0);
-    addToRunTimeSelectionTable
-    (
-        wallHeatFluxModel,
-        wall,
-        dictionary
-    );
-}
-}
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-void Foam::wallHeatFluxModels::wall::writeFileHeader(Ostream& os)
-{
-    writeHeader(os, "Wall heat-flux");
-    writeCommented(os, "Time");
-    writeTabbed(os, "patch");
-    writeTabbed(os, "min");
-    writeTabbed(os, "max");
-    writeTabbed(os, "integral");
-    os  << endl;
-
-    writtenHeader_ = true;
-}
-
-
-void Foam::wallHeatFluxModels::wall::calcHeatFlux
-(
-    const volScalarField& alpha,
-    const volScalarField& he,
-    volScalarField& wallHeatFlux
-)
-{
-    volScalarField::Boundary& wallHeatFluxBf = wallHeatFlux.boundaryFieldRef();
-
-    const volScalarField::Boundary& heBf = he.boundaryField();
-
-    const volScalarField::Boundary& alphaBf = alpha.boundaryField();
-
-    const labelHashSet& patches = patchSet();
-    for (const label patchi : patches)
-    {
-        wallHeatFluxBf[patchi] = alphaBf[patchi]*heBf[patchi].snGrad();
-    }
-
-
-    const auto* qrPtr = mesh().cfindObject<volScalarField>(qrName());
-
-    if (qrPtr)
-    {
-        const volScalarField::Boundary& radHeatFluxBf = qrPtr->boundaryField();
-
-        for (const label patchi : patches)
-        {
-            wallHeatFluxBf[patchi] -= radHeatFluxBf[patchi];
-        }
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::wallHeatFluxModels::wall::wall
-(
-    const dictionary& dict,
-    const fvMesh& mesh,
-    const word& name,
-    const word objName,
-    functionObjects::stateFunctionObject& state
-)
-:
-    wallHeatFluxModel(dict, mesh, name, objName, state)
-{
-    auto* wallHeatFluxPtr
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                objName,
-                mesh.time().timeName(),
-                mesh
-            ),
-            mesh,
-            dimensionedScalar(dimMass/pow3(dimTime), Zero)
-        )
-    );
-
-    mesh.objectRegistry::store(wallHeatFluxPtr);
-}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-bool Foam::wallHeatFluxModels::wall::read(const dictionary& dict)
-{
-    if (!wallHeatFluxModel::read(dict))
-    {
-        return false;
-    }
-
-
-    qrName_ = dict.getOrDefault<word>("qr", "qr");
-
-
-    Info<< state().type() << " " << state().name() << ":" << nl;
-
-    patchSet_ = mesh().boundaryMesh().patchSet
-    (
-        dict.getOrDefault<wordRes>("patches", wordRes())
-    );
-
-    const polyBoundaryMesh& pbm = mesh().boundaryMesh();
-    if (patchSet_.empty())
-    {
-        forAll(pbm, patchi)
-        {
-            if (isA<wallPolyPatch>(pbm[patchi]))
-            {
-                patchSet_.insert(patchi);
-            }
-        }
-
-        Info<< "    processing all wall patches" << nl << endl;
-    }
-    else
-    {
-        Info<< "    processing wall patches: " << nl;
-        labelHashSet filteredPatchSet;
-        for (const label patchi : patchSet_)
-        {
-            if (isA<wallPolyPatch>(pbm[patchi]))
-            {
-                filteredPatchSet.insert(patchi);
-                Info<< "        " << pbm[patchi].name() << endl;
-            }
-            else
-            {
-                WarningInFunction
-                    << "Requested wall heat-flux on non-wall boundary "
-                    << "type patch: " << pbm[patchi].name() << endl;
-            }
-        }
-
-        Info<< endl;
-
-        patchSet_ = filteredPatchSet;
-    }
-
-
-    if (writeFile::canResetFile())
-    {
-        writeFile::resetFile(objName());
-    }
-
-    if (writeFile::canWriteHeader())
-    {
-        writeFileHeader(file());
-    }
-
-
-    return true;
-}
-
-
-bool Foam::wallHeatFluxModels::wall::execute()
-{
-    auto& wallHeatFlux = mesh().lookupObjectRef<volScalarField>(objName());
-
-    if
-    (
-        mesh().foundObject<compressible::turbulenceModel>
-        (
-            turbulenceModel::propertiesName
-        )
-    )
-    {
-        const compressible::turbulenceModel& turbModel =
-            mesh().lookupObject<compressible::turbulenceModel>
-            (
-                turbulenceModel::propertiesName
-            );
-
-        calcHeatFlux
-        (
-            turbModel.alphaEff()(),
-            turbModel.transport().he(),
-            wallHeatFlux
-        );
-    }
-    else if (mesh().foundObject<fluidThermo>(fluidThermo::dictName))
-    {
-        const fluidThermo& thermo =
-            mesh().lookupObject<fluidThermo>(fluidThermo::dictName);
-
-        calcHeatFlux
-        (
-            thermo.alpha(),
-            thermo.he(),
-            wallHeatFlux
-        );
-    }
-    else if (mesh().foundObject<solidThermo>(solidThermo::dictName))
-    {
-        const solidThermo& thermo =
-            mesh().lookupObject<solidThermo>(solidThermo::dictName);
-
-        calcHeatFlux(thermo.alpha(), thermo.he(), wallHeatFlux);
-    }
-    else if
-    (
-        mesh().foundObject<multiphaseInterSystem>
-            (multiphaseInterSystem::phasePropertiesName)
-    )
-    {
-        const auto& thermo = mesh().lookupObject<multiphaseInterSystem>
-        (
-            multiphaseInterSystem::phasePropertiesName
-        );
-
-        calcHeatFlux(thermo.kappaEff()(), thermo.T(), wallHeatFlux);
-    }
-    else
-    {
-        FatalErrorInFunction
-            << "Unable to find compressible turbulence model in the "
-            << "database" << exit(FatalError);
-    }
-
-    const fvPatchList& patches = mesh().boundary();
-
-    const surfaceScalarField::Boundary& magSf = mesh().magSf().boundaryField();
-
-    const labelHashSet& patchset = patchSet();
-    for (const label patchi : patchset)
-    {
-        const fvPatch& pp = patches[patchi];
-
-        const scalarField& hfp = wallHeatFlux.boundaryField()[patchi];
-
-        const scalar minHfp = gMin(hfp);
-        const scalar maxHfp = gMax(hfp);
-        const scalar integralHfp = gSum(magSf[patchi]*hfp);
-
-        if (Pstream::master())
-        {
-            writeCurrentTime(file());
-
-            file()
-                << token::TAB << pp.name()
-                << token::TAB << minHfp
-                << token::TAB << maxHfp
-                << token::TAB << integralHfp
-                << endl;
-        }
-
-        if (state().log)
-        {
-            Info<< "    min/max/integ(" << pp.name() << ") = "
-                << minHfp << ", " << maxHfp << ", " << integralHfp << endl;
-        }
-
-        state().setResult("min(" + pp.name() + ")", minHfp);
-        state().setResult("max(" + pp.name() + ")", maxHfp);
-        state().setResult("int(" + pp.name() + ")", integralHfp);
-    }
-
-
-    return true;
-}
-
-
-bool Foam::wallHeatFluxModels::wall::write()
-{
-    const auto& wallHeatFlux =
-        mesh().lookupObject<volScalarField>(objName());
-
-    if (state().log)
-    {
-        Info<< "    writing field " << wallHeatFlux.name() << endl;
-    }
-
-    wallHeatFlux.write();
-
-    return true;
-}
-
-
-// ************************************************************************* //

src/functionObjects/field/wallHeatFlux/wallHeatFluxModels/wall/wall.H

@@ -1,176 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2025 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::wallHeatFluxModels::wall
-
-Description
-    This model computes the wall-heat flux at selected wall patches.
-
-Usage
-    Minimal example by using \c system/controlDict.functions:
-    \verbatim
-    wallHeatFlux1
-    {
-        // Mandatory entries
-        type        wallHeatFlux;
-        libs        (fieldFunctionObjects);
-
-        model       wall;
-
-        // Optional entries
-        patches     (<patch1> ... <patchN>); // (wall1 "(wall2|wall3)");
-        qr          <word>;
-
-        // Inherited entries
-        ...
-    }
-    \endverbatim
-
-    where the entries mean:
-    \table
-      Property     | Description                           | Type | Reqd | Deflt
-      type         | Type name: wallHeatFlux               | word |  yes  | -
-      libs         | Library name: fieldFunctionObjects    | word |  yes  | -
-      patches    | Names of operand patches   | wordList | no | all wall patches
-      qr         | Name of radiative heat flux field       | word | no | qr
-    \endtable
-
-    The inherited entries are elaborated in:
-     - \link fvMeshFunctionObject.H \endlink
-     - \link writeFile.H \endlink
-
-Note
-  - The model \c wall corresponds to \c wallHeatFlux in OpenFOAM v2506 and
-    earlier versions.
-
-SourceFiles
-    wall.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef wallHeatFluxModels_wall_H
-#define wallHeatFluxModels_wall_H
-
-#include "wallHeatFluxModel.H"
-#include "writeFile.H"
-#include "volFieldsFwd.H"
-#include "HashSet.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace wallHeatFluxModels
-{
-
-/*---------------------------------------------------------------------------*\
-                   Class wall Declaration
-\*---------------------------------------------------------------------------*/
-
-class wall
-:
-    public wallHeatFluxModel
-{
-    // Private Data
-
-        //- Name of the radiative heat-flux field
-        word qrName_;
-
-        //- List of names of wall patches to process
-        labelHashSet patchSet_;
-
-
-    // Private Member Functions
-
-        //- Write file-header information
-        void writeFileHeader(Ostream& os);
-
-        //- Calculate the wall heat-flux
-        void calcHeatFlux
-        (
-            const volScalarField& alpha,
-            const volScalarField& he,
-            volScalarField& wallHeatFlux
-        );
-
-
-public:
-
-    //- Runtime type information
-    TypeName("wall");
-
-
-    // Constructors
-
-        //- Construct from components
-        wall
-        (
-            const dictionary& dict,
-            const fvMesh& mesh,
-            const word& name,
-            const word objName,
-            functionObjects::stateFunctionObject& state
-        );
-
-
-    //- Destructor
-    virtual ~wall() = default;
-
-
-    // Member Functions
-
-    // Access
-
-        //- Return const reference to name of radiative heat-flux field
-        const word& qrName() const noexcept { return qrName_; }
-
-        //- Return const reference to patches to process
-        const labelHashSet& patchSet() const noexcept { return patchSet_; }
-
-
-    // Evaluation
-
-        //- Read the settings
-        virtual bool read(const dictionary& dict);
-
-        //- Calculate the heat-flux data
-        virtual bool execute();
-
-        //- Write the heat-flux data
-        virtual bool write();
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace wallHeatFluxModels
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/functionObjects/field/wallHeatFlux/wallHeatFluxModels/wallHeatFluxModel/wallHeatFluxModel.C

@@ -1,73 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2025 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 "wallHeatFluxModel.H"
-#include "fvMesh.H"
-#include "wallPolyPatch.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTypeNameAndDebug(wallHeatFluxModel, 0);
-    defineRunTimeSelectionTable(wallHeatFluxModel, dictionary);
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::wallHeatFluxModel::wallHeatFluxModel
-(
-    const dictionary& dict,
-    const fvMesh& mesh,
-    const word& name,
-    const word objName,
-    functionObjects::stateFunctionObject& state
-)
-:
-    functionObjects::writeFile(mesh, name, objName, dict),
-    mesh_(mesh),
-    state_(state),
-    objName_(objName)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-bool Foam::wallHeatFluxModel::read(const dictionary& dict)
-{
-    if (!writeFile::read(dict))
-    {
-        return false;
-    }
-
-    return true;
-}
-
-
-// ************************************************************************* //

src/functionObjects/field/wallHeatFlux/wallHeatFluxModels/wallHeatFluxModel/wallHeatFluxModel.H

@@ -1,181 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2025 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/>.
-
-Namespace
-    Foam::wallHeatFluxModels
-
-Description
-    A namespace for various heat-flux model implementations.
-
-Class
-    Foam::wallHeatFluxModel
-
-Description
-    A base class for heat-flux models.
-
-SourceFiles
-    wallHeatFluxModel.C
-    wallHeatFluxModelNew.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef Foam_wallHeatFluxModel_H
-#define Foam_wallHeatFluxModel_H
-
-#include "writeFile.H"
-#include "volFieldsFwd.H"
-#include "HashSet.H"
-#include "stateFunctionObject.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward Declarations
-class fvMesh;
-
-/*---------------------------------------------------------------------------*\
-                   Class wallHeatFluxModel Declaration
-\*---------------------------------------------------------------------------*/
-
-class wallHeatFluxModel
-:
-    public functionObjects::writeFile
-{
-    // Private Data
-
-        //- Const reference to the mesh
-        const fvMesh& mesh_;
-
-        //- Reference to the state function object
-        functionObjects::stateFunctionObject& state_;
-
-        //- Name of the function object
-        const word objName_;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("wallHeatFluxModel");
-
-
-    // Declare runtime constructor selection table
-
-        declareRunTimeSelectionTable
-        (
-            autoPtr,
-            wallHeatFluxModel,
-            dictionary,
-            (
-                const dictionary& dict,
-                const fvMesh& mesh,
-                const word& name,
-                const word objName,
-                functionObjects::stateFunctionObject& state
-            ),
-            (dict, mesh, name, objName, state)
-        );
-
-
-    // Selectors
-
-        //- Return a reference to the selected heat-flux model
-        static autoPtr<wallHeatFluxModel> New
-        (
-            const dictionary& dict,
-            const fvMesh& mesh,
-            const word& name,
-            const word objName,
-            functionObjects::stateFunctionObject& state
-        );
-
-
-    // Generated Methods
-
-        //- No copy construct
-        wallHeatFluxModel(const wallHeatFluxModel&) = delete;
-
-        //- No copy assignment
-        void operator=(const wallHeatFluxModel&) = delete;
-
-
-    // Constructors
-
-        //- Construct from components
-        wallHeatFluxModel
-        (
-            const dictionary& dict,
-            const fvMesh& mesh,
-            const word& name,
-            const word objName,
-            functionObjects::stateFunctionObject& state
-        );
-
-
-    //- Destructor
-    virtual ~wallHeatFluxModel() = default;
-
-
-    // Member Functions
-
-    // Access
-
-        //- Return const reference to the mesh
-        const fvMesh& mesh() const noexcept { return mesh_; }
-
-        //- Return const reference to the state function object
-        functionObjects::stateFunctionObject& state() const noexcept
-        {
-            return state_;
-        }
-
-        //- Return const reference to the function-object name
-        const word& objName() const noexcept { return objName_; }
-
-
-    // Evaluation
-
-        //- Read the settings
-        virtual bool read(const dictionary& dict) = 0;
-
-        //- Calculate the heat-flux data
-        virtual bool execute() = 0;
-
-        //- Write the heat-flux data
-        virtual bool write() = 0;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/functionObjects/field/wallHeatFlux/wallHeatFluxModels/wallHeatFluxModel/wallHeatFluxModelNew.C

@@ -1,67 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2025 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 "wallHeatFluxModel.H"
-#include "fvMesh.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-Foam::autoPtr<Foam::wallHeatFluxModel> Foam::wallHeatFluxModel::New
-(
-    const dictionary& dict,
-    const fvMesh& mesh,
-    const word& name,
-    const word objName,
-    functionObjects::stateFunctionObject& state
-)
-{
-    const word modelType(dict.getOrDefault<word>("model", "wall"));
-
-    Info<< "Selecting heat-flux model: " << modelType << endl;
-
-    auto* ctorPtr = dictionaryConstructorTable(modelType);
-
-    if (!ctorPtr)
-    {
-        FatalIOErrorInLookup
-        (
-            dict,
-            "wallHeatFluxModel",
-            modelType,
-            *dictionaryConstructorTablePtr_
-        ) << exit(FatalIOError);
-    }
-
-    return autoPtr<wallHeatFluxModel>
-    (
-        ctorPtr(dict, mesh, name, objName, state)
-    );
-}
-
-
-// ************************************************************************* //

src/fused/finiteVolume/fusedGaussGrad.C

@@ -226,6 +226,94 @@ Foam::fv::fusedGaussGrad<Type>::calcGrad
 }
 
 
+template<class Type>
+void Foam::fv::fusedGaussGrad<Type>::calcGrad
+(
+    GeometricField
+    <
+        typename outerProduct<vector, Type>::type,
+        fvPatchField,
+        volMesh
+    >& gGrad,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+) const
+{
+    typedef typename outerProduct<vector, Type>::type GradType;
+
+    const fvMesh& mesh = vf.mesh();
+
+    DebugPout<< "fusedGaussGrad<Type>::calcGrad on " << vf.name()
+        << " into " << gGrad.name() << endl;
+
+    gGrad = dimensioned<GradType>(vf.dimensions()/dimLength, Zero);
+
+    if (this->tinterpScheme_().corrected())
+    {
+        const auto tfaceCorr(this->tinterpScheme_().correction(vf));
+        auto& faceCorr = tfaceCorr();
+
+        DebugPout<< "fusedGaussGrad<Type>::calcGrad corrected interpScheme "
+            << this->tinterpScheme_().type() << endl;
+
+        const auto interpolate = [&]
+        (
+            const vector& area,
+            const scalar lambda,
+
+            const Type& ownVal,
+            const Type& neiVal,
+
+            const Type& correction
+
+        ) -> GradType
+        {
+            return area*((lambda*(ownVal - neiVal) + neiVal) + correction);
+        };
+
+        fvc::surfaceSum
+        (
+            this->tinterpScheme_().weights(vf),
+            vf,
+            faceCorr,
+            interpolate,
+            gGrad,
+            false
+        );
+    }
+    else
+    {
+        DebugPout<< "fusedGaussGrad<Type>::calcGrad uncorrected interpScheme "
+            << this->tinterpScheme_().type() << endl;
+
+        const auto interpolate = [&]
+        (
+            const vector& area,
+            const scalar lambda,
+            const Type& ownVal,
+            const Type& neiVal
+        ) -> GradType
+        {
+            return area*(lambda*(ownVal - neiVal) + neiVal);
+        };
+
+        fvc::surfaceSum
+        (
+            tinterpScheme_().weights(vf),
+            vf,
+            interpolate,
+            gGrad,
+            false
+        );
+    }
+
+    gGrad.primitiveFieldRef() /= mesh.V();
+
+    gGrad.correctBoundaryConditions();
+
+    correctBoundaryConditions(vf, gGrad);
+}
+
+
 template<class Type>
 template<class GradType>
 void Foam::fv::fusedGaussGrad<Type>::correctBoundaryConditions

src/fused/finiteVolume/fusedGaussGrad.H

@@ -145,6 +145,14 @@ public:
             const word& name
         ) const;
 
+        //- Calculate the grad of the given field into supplied field
+        virtual void calcGrad
+        (
+            GeometricField
+            <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>& res,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const;
+
         //- Correct the boundary values of the gradient using the patchField
         //- snGrad functions
         template<class GradType>

src/fused/finiteVolume/fusedGaussLaplacianScheme.C

@@ -45,6 +45,87 @@ namespace fv
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+template<class Type, class GType>
+template<class E1, class E2>
+void fusedGaussLaplacianScheme<Type, GType>::fvmCorrection
+(
+    fvMatrix<Type>& fvm,
+    const dimensionSet& gammaDim,
+    const E1& gammaMagSf,
+    const E2& corr
+)
+{
+    typedef GeometricField<Type, fvsPatchField, surfaceMesh> surfaceType;
+
+    const auto& vf = fvm.psi();
+    const fvMesh& mesh = vf.mesh();
+    const auto V = mesh.V().expr();
+
+    if (mesh.fluxRequired(vf.name()))
+    {
+        fvm.faceFluxCorrectionPtr() = std::make_unique<surfaceType>
+        (
+            IOobject
+            (
+                "faceFluxCorr",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                IOobject::NO_REGISTER
+            ),
+            mesh,
+            gammaDim
+           *mesh.magSf().dimensions()
+           *corr.data().dimensions()
+        );
+        auto& faceFluxCorr = *fvm.faceFluxCorrectionPtr();
+
+        faceFluxCorr = corr*gammaMagSf;
+
+        fvm.source() =
+            fvm.source().expr()
+          - (
+                fvc::div
+                (
+                    faceFluxCorr
+                )().primitiveField().expr()
+               *V
+            );
+    }
+    else
+    {
+        surfaceType faceFluxCorr
+        (
+            IOobject
+            (
+                "faceFluxCorr",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                IOobject::NO_REGISTER
+            ),
+            mesh,
+            gammaDim
+           *mesh.magSf().dimensions()
+           *corr.data().dimensions()
+        );
+        faceFluxCorr = corr*gammaMagSf;
+
+        fvm.source() =
+            fvm.source().expr()
+          - (
+                fvc::div
+                (
+                    faceFluxCorr
+                )().primitiveField().expr()
+               *V
+            );
+    }
+}
+
+
 template<class Type, class GType>
 tmp<fvMatrix<Type>>
 fusedGaussLaplacianScheme<Type, GType>::fvmLaplacianUncorrected

src/fused/finiteVolume/fusedGaussLaplacianScheme.H

@@ -65,6 +65,17 @@ class fusedGaussLaplacianScheme
 {
     // Private Member Functions
 
+        //- Helper for calculate source correction. Partially expression
+        //  templated. Move to fvMatrixExpression.H?
+        template<class E1, class E2>
+        static void fvmCorrection
+        (
+            fvMatrix<Type>& fvm,
+            const dimensionSet& gammaDim,
+            const E1& gammaMagSf,
+            const E2& corr
+        );
+
         //- See gaussLaplacianScheme
         tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> gammaSnGradCorr
         (
@@ -198,6 +209,8 @@ defineFvmLaplacianScalarGamma(sphericalTensor);
 defineFvmLaplacianScalarGamma(symmTensor);
 defineFvmLaplacianScalarGamma(tensor);
 
+#undef defineFvmLaplacianScalarGamma
+
 // Unweighted laplacian
 template<>
 tmp<GeometricField<scalar, fvPatchField, volMesh>>
@@ -227,20 +240,6 @@ fusedGaussLaplacianScheme<vector, scalar>::fvcLaplacian
     const GeometricField<scalar, fvPatchField, volMesh>&,
     const GeometricField<vector, fvPatchField, volMesh>&
 );
-template<>
-tmp<fvMatrix<scalar>>
-fusedGaussLaplacianScheme<scalar, scalar>::fvmLaplacian
-(
-    const GeometricField<scalar, fvPatchField, volMesh>&,
-    const GeometricField<scalar, fvPatchField, volMesh>&
-);
-template<>
-tmp<fvMatrix<vector>>
-fusedGaussLaplacianScheme<vector, scalar>::fvmLaplacian
-(
-    const GeometricField<scalar, fvPatchField, volMesh>&,
-    const GeometricField<vector, fvPatchField, volMesh>&
-);
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/fused/finiteVolume/fusedGaussLaplacianSchemes.C

@@ -28,6 +28,7 @@ License
 
 #include "fusedGaussLaplacianScheme.H"
 #include "fvMesh.H"
+#include "fvMatrixExpression.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -40,61 +41,81 @@ Foam::tmp<Foam::fvMatrix<Foam::Type>>                                          \
 Foam::fv::fusedGaussLaplacianScheme<Foam::Type, Foam::scalar>::                \
 fvmLaplacian                                                                   \
 (                                                                              \
-    const GeometricField<scalar, fvsPatchField, surfaceMesh>& gamma,           \
+    const GeometricField<scalar, fvPatchField, volMesh>& gamma,                \
     const GeometricField<Type, fvPatchField, volMesh>& vf                      \
 )                                                                              \
 {                                                                              \
     DebugPout<< "fusedGaussLaplacianScheme::fvmLaplacian on " << vf.name()     \
-        << " with scalar gamma " << gamma.name() << endl;                      \
+        << " with scalar gamma " << gamma.name() << " and ET" << endl;         \
                                                                                \
     const fvMesh& mesh = this->mesh();                                         \
                                                                                \
-    GeometricField<scalar, fvsPatchField, surfaceMesh> gammaMagSf              \
-    (                                                                          \
-        gamma*mesh.magSf()                                                     \
-    );                                                                         \
+    const auto weights = this->tinterpGammaScheme_().weights(gamma).expr();    \
+    const auto gammaMagSf =                                                    \
+        Expression::interpolate(gamma.expr(), weights, mesh)                   \
+      * mesh.magSf().expr();                                                   \
+    /* For compatibility with linearInterpolate : avoid orientation. TBD. */   \
+    const_cast<orientedType&>(gammaMagSf.oriented()).setOriented(false);       \
+    const auto deltaCoeffs = this->tsnGradScheme_().deltaCoeffs(vf).expr();    \
                                                                                \
-    tmp<fvMatrix<Type>> tfvm = fvmLaplacianUncorrected                         \
+    tmp<fvMatrix<Type>> tfvm                                                   \
     (                                                                          \
-        gammaMagSf,                                                            \
-        this->tsnGradScheme_().deltaCoeffs(vf),                                \
-        vf                                                                     \
+        new fvMatrix<Type>                                                     \
+        (                                                                      \
+            vf,                                                                \
+            gamma.dimensions()*mesh.magSf().dimensions()*vf.dimensions()       \
+        )                                                                      \
     );                                                                         \
     fvMatrix<Type>& fvm = tfvm.ref();                                          \
                                                                                \
+    Expression::fvmLaplacianUncorrected(fvm, gammaMagSf, deltaCoeffs);         \
+                                                                               \
     if (this->tsnGradScheme_().corrected())                                    \
     {                                                                          \
-        if (mesh.fluxRequired(vf.name()))                                      \
-        {                                                                      \
-            fvm.faceFluxCorrectionPtr() = std::make_unique                     \
-            <                                                                  \
-                GeometricField<Type, fvsPatchField, surfaceMesh>               \
-            >                                                                  \
-            (                                                                  \
-                gammaMagSf*this->tsnGradScheme_().correction(vf)               \
-            );                                                                 \
-                                                                               \
-            fvm.source() -=                                                    \
-                mesh.V()*                                                      \
-                fvc::div                                                       \
-                (                                                              \
-                    *fvm.faceFluxCorrectionPtr()                               \
-                )().primitiveField();                                          \
-        }                                                                      \
-        else                                                                   \
-        {                                                                      \
-            fvm.source() -=                                                    \
-                mesh.V()*                                                      \
-                fvc::div                                                       \
-                (                                                              \
-                    gammaMagSf*this->tsnGradScheme_().correction(vf)           \
-                )().primitiveField();                                          \
-        }                                                                      \
+        const auto corr(this->tsnGradScheme_().correction(vf).expr());         \
+        fvmCorrection(fvm, gamma.dimensions(), gammaMagSf, corr);              \
     }                                                                          \
                                                                                \
     return tfvm;                                                               \
 }                                                                              \
                                                                                \
+template<>                                                                     \
+Foam::tmp<Foam::fvMatrix<Foam::Type>>                                          \
+Foam::fv::fusedGaussLaplacianScheme<Foam::Type, Foam::scalar>::                \
+fvmLaplacian                                                                   \
+(                                                                              \
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>& gamma,           \
+    const GeometricField<Type, fvPatchField, volMesh>& vf                      \
+)                                                                              \
+{                                                                              \
+    DebugPout<< "fusedGaussLaplacianScheme::fvmLaplacian on " << vf.name()     \
+        << " with interpolated gamma " << gamma.name() << " and ET" << endl;   \
+                                                                               \
+    const fvMesh& mesh = this->mesh();                                         \
+                                                                               \
+    const auto gammaMagSf = gamma.expr()* mesh.magSf().expr();                 \
+    const auto deltaCoeffs = this->tsnGradScheme_().deltaCoeffs(vf).expr();    \
+                                                                               \
+    tmp<fvMatrix<Type>> tfvm                                                   \
+    (                                                                          \
+        new fvMatrix<Type>                                                     \
+        (                                                                      \
+            vf,                                                                \
+            gamma.dimensions()*mesh.magSf().dimensions()*vf.dimensions()       \
+        )                                                                      \
+    );                                                                         \
+    fvMatrix<Type>& fvm = tfvm.ref();                                          \
+                                                                               \
+    Expression::fvmLaplacianUncorrected(fvm, gammaMagSf, deltaCoeffs);         \
+                                                                               \
+    if (this->tsnGradScheme_().corrected())                                    \
+    {                                                                          \
+        const auto corr(this->tsnGradScheme_().correction(vf).expr());         \
+        fvmCorrection(fvm, gamma.dimensions(), gammaMagSf, corr);              \
+    }                                                                          \
+                                                                               \
+    return tfvm;                                                               \
+}                                                                              \
                                                                                \
 template<>                                                                     \
 Foam::tmp<Foam::GeometricField<Foam::Type, Foam::fvPatchField, Foam::volMesh>> \
@@ -500,38 +521,6 @@ Foam::fv::fusedGaussLaplacianScheme<Foam::vector, Foam::scalar>::fvcLaplacian
 }
 
 
-template<>
-Foam::tmp<Foam::fvMatrix<Foam::scalar>>
-Foam::fv::fusedGaussLaplacianScheme<Foam::scalar, Foam::scalar>::fvmLaplacian
-(
-    const GeometricField<scalar, fvPatchField, volMesh>& gamma,
-    const GeometricField<scalar, fvPatchField, volMesh>& vf
-)
-{
-    // TBD
-    DebugPout
-        << "fusedGaussLaplacianScheme<scalar, scalar>::fvmLaplacian"
-        << " on " << vf.name() << " with gamma " << gamma.name() << endl;
-    return fvmLaplacian(this->tinterpGammaScheme_().interpolate(gamma)(), vf);
-}
-
-
-template<>
-Foam::tmp<Foam::fvMatrix<Foam::vector>>
-Foam::fv::fusedGaussLaplacianScheme<Foam::vector, Foam::scalar>::fvmLaplacian
-(
-    const GeometricField<scalar, fvPatchField, volMesh>& gamma,
-    const GeometricField<vector, fvPatchField, volMesh>& vf
-)
-{
-    // TBD
-    DebugPout
-        << "fusedGaussLaplacianScheme<vector, scalar>::fvmLaplacian"
-        << " on " << vf.name() << " with gamma " << gamma.name() << endl;
-    return fvmLaplacian(this->tinterpGammaScheme_().interpolate(gamma)(), vf);
-}
-
-
 template<>
 Foam::tmp
 <

src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/dynamicContactAngleForce/dynamicContactAngleForce.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020-2023 OpenCFD Ltd.
+    Copyright (C) 2020-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -26,9 +26,9 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "dynamicContactAngleForce.H"
-#include "addToRunTimeSelectionTable.H"
 #include "Function1.H"
 #include "distributionModel.H"
+#include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -59,6 +59,7 @@ dynamicContactAngleForce::dynamicContactAngleForce
 )
 :
     contactAngleForce(typeName, film, dict),
+    thetaPtr_(nullptr),
     U_vs_thetaPtr_
     (
         Function1<scalar>::NewIfPresent
@@ -80,46 +81,56 @@ dynamicContactAngleForce::dynamicContactAngleForce
         )
     ),
     rndGen_(label(0)),
-    distribution_
-    (
-        distributionModel::New
-        (
-            coeffDict_.subDict("distribution"),
-            rndGen_
-        )
-    )
+    distributionPtr_(nullptr)
 {
     if (U_vs_thetaPtr_ && T_vs_thetaPtr_)
     {
         FatalIOErrorInFunction(dict)
-            << "Entries Utheta and Ttheta could not be used together"
+            << "Only one of Utheta or Ttheta should be provided; "
+            << "both inputs cannot be used together."
             << abort(FatalIOError);
     }
+
+
+    thetaPtr_.emplace
+    (
+        IOobject
+        (
+            IOobject::scopedName(typeName, "theta"),
+            film.regionMesh().time().timeName(),
+            film.regionMesh().thisDb(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        film.regionMesh(),
+        dimensionedScalar(dimless, Zero)
+    );
+
+
+    if (coeffDict_.findEntry("distribution"))
+    {
+        distributionPtr_ = distributionModel::New
+        (
+            coeffDict_.subDict("distribution"),
+            rndGen_
+        );
+    }
 }
 
 
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+dynamicContactAngleForce::~dynamicContactAngleForce()
+{}  // distributionModel was forward declared
+
+
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 tmp<areaScalarField> dynamicContactAngleForce::theta() const
 {
-    auto ttheta = tmp<areaScalarField>::New
-    (
-        IOobject
-        (
-            IOobject::scopedName(typeName, "theta"),
-            film().regionMesh().time().timeName(),
-            film().regionMesh().thisDb(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            IOobject::NO_REGISTER
-        ),
-        film().regionMesh(),
-        dimensionedScalar(dimless, Zero)
-    );
-    areaScalarField& theta = ttheta.ref();
+    areaScalarField& theta = thetaPtr_.ref();
     scalarField& thetai = theta.ref();
 
-
     if (U_vs_thetaPtr_)
     {
         // Initialize with the function of film speed
@@ -136,13 +147,16 @@ tmp<areaScalarField> dynamicContactAngleForce::theta() const
         thetai = T_vs_thetaPtr_->value(T());
     }
 
-    // Add the stochastic perturbation
-    forAll(thetai, facei)
+    if (distributionPtr_)
     {
-        thetai[facei] += distribution_->sample();
+        // Add the stochastic perturbation
+        forAll(thetai, facei)
+        {
+            thetai[facei] += distributionPtr_->sample();
+        }
     }
 
-    return ttheta;
+    return thetaPtr_();
 }
 
 

src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/dynamicContactAngleForce/dynamicContactAngleForce.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020-2022 OpenCFD Ltd.
+    Copyright (C) 2020-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -112,6 +112,9 @@ class dynamicContactAngleForce
 {
     // Private Data
 
+        //- Contact-angle field in degrees
+        mutable autoPtr<areaScalarField> thetaPtr_;
+
         //- Contact angle as a function of film speed
         autoPtr<Function1<scalar>> U_vs_thetaPtr_;
 
@@ -121,17 +124,8 @@ class dynamicContactAngleForce
         //- Random number generator
         Random rndGen_;
 
-        //- Parcel size PDF model
-        const autoPtr<distributionModel> distribution_;
-
-
-    // Private Member Functions
-
-        //- No copy construct
-        dynamicContactAngleForce(const dynamicContactAngleForce&) = delete;
-
-        //- No copy assignment
-        void operator=(const dynamicContactAngleForce&) = delete;
+        //- Stochastic perturbation model for contact angle
+        autoPtr<distributionModel> distributionPtr_;
 
 
 protected:
@@ -148,7 +142,7 @@ public:
 
     // Constructors
 
-        //- Construct from surface film model
+        //- Construct from surface film model and dictionary
         dynamicContactAngleForce
         (
             liquidFilmBase& film,
@@ -157,7 +151,7 @@ public:
 
 
     //- Destructor
-    virtual ~dynamicContactAngleForce() = default;
+    virtual ~dynamicContactAngleForce();
 };
 
 

tutorials/incompressible/simpleFoam/pitzDaily_fused/Allrun

@@ -19,4 +19,12 @@ sed "s/GAUSS/fusedGauss/g" system/fvSchemes.template > system/fvSchemes
 runParallel -s fusedGauss $(getApplication)
 
 
+echo "Updating fvSolution to cache gradients"
+sed 's/GRADIENT/"grad(.*)"/g' system/fvSolution.template > system/fvSolution
+
+runParallel -s fusedGaussCached $(getApplication)
+
+# Restore
+cp system/fvSolution.template system/fvSolution
+
 #------------------------------------------------------------------------------

tutorials/incompressible/simpleFoam/pitzDaily_fused/system/controlDict

@@ -16,6 +16,15 @@ FoamFile
 
 libs            (fusedFiniteVolume);
 
+/*
+DebugSwitches
+{
+    fusedGauss  1;
+    Gauss       1;
+    solution    1;
+}
+*/
+
 application     simpleFoam;
 
 startFrom       startTime;

tutorials/incompressible/simpleFoam/pitzDaily_fused/system/fvSolution

@@ -55,5 +55,11 @@ relaxationFactors
     }
 }
 
+cache
+{
+    // Cache all grads
+    GRADIENT;
+}
+
 
 // ************************************************************************* //

tutorials/incompressible/simpleFoam/pitzDaily_fused/system/fvSolution.template

@@ -0,0 +1,65 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2506                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p
+    {
+        solver          GAMG;
+        tolerance       1e-06;
+        relTol          0.1;
+        smoother        GaussSeidel;
+    }
+
+    "(U|k|epsilon|omega|f|v2)"
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-05;
+        relTol          0.1;
+    }
+}
+
+SIMPLE
+{
+    nNonOrthogonalCorrectors 0;
+    consistent      yes;
+
+    residualControl
+    {
+        p               1e-2;
+        U               1e-3;
+        "(k|epsilon|omega|f|v2)" 1e-3;
+    }
+}
+
+relaxationFactors
+{
+    equations
+    {
+        U               0.9; // 0.9 is more stable but 0.95 more convergent
+        ".*"            0.9; // 0.9 is more stable but 0.95 more convergent
+    }
+}
+
+cache
+{
+    // Cache all grads
+    GRADIENT;
+}
+
+
+// ************************************************************************* //