ENH: dynamicContactAngle: enable zonal input using persistent field storage

Temporary field creation is replaced with persistent field storage,
so that contact-angle can be input per zone using setFields utility.

applications/test/expressionTemplates/Make/files

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

applications/test/expressionTemplates/Make/options

@@ -1,11 +0,0 @@
-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

@@ -1,126 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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

@@ -1,451 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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,6 +55,8 @@ 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,7 +45,6 @@ SourceFiles
 
 #include "autoPtr.H"
 #include "UList.H"
-#include "SLListFwd.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -66,6 +65,7 @@ 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,37 +395,6 @@ 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,11 +80,6 @@ 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
@@ -702,39 +697,6 @@ 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,12 +60,6 @@ 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
 \*---------------------------------------------------------------------------*/
@@ -455,28 +449,6 @@ 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

@@ -1,252 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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

@@ -1,145 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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

@@ -1,148 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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,37 +582,6 @@ 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,124 +1675,6 @@ 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,7 +42,6 @@ SourceFiles
 #define Foam_GeometricField_H
 
 #include "GeometricBoundaryField.H"
-#include "GeometricFieldExpression.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -1143,79 +1142,6 @@ 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,17 +59,8 @@ 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
@@ -361,33 +352,6 @@ 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/Pstream/mpi/UPstreamRequest.C

@@ -855,7 +855,7 @@ bool Foam::UPstream::finishedRequest(const label i)
     // This allows MPI to progress behind the scenes if it wishes.
 
     int flag = 0;
-    if (i >= 0 && i < PstreamGlobals::outstandingRequests_.size())
+    if (i < 0 || i >= PstreamGlobals::outstandingRequests_.size())
     {
         auto& request = PstreamGlobals::outstandingRequests_[i];
 

src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.C

@@ -218,11 +218,9 @@ tmp<vectorField> atmBoundaryLayer::U(const vectorField& pCf) const
     const scalar groundMin = zDir() & ppMin_;
 
     // (YGCJ:Table 1, RH:Eq. 6, HW:Eq. 5)
-    scalarField zEff(max((zDir() & pCf) - groundMin - d + z0, z0));
-
     scalarField Un
     (
-        (Ustar(z0)/kappa_)*log(zEff/z0)
+        (Ustar(z0)/kappa_)*log(((zDir() & pCf) - groundMin - d + z0)/z0)
     );
 
     return flowDir()*Un;
@@ -237,9 +235,9 @@ tmp<scalarField> atmBoundaryLayer::k(const vectorField& pCf) const
     const scalar groundMin = zDir() & ppMin_;
 
     // (YGCJ:Eq. 21; RH:Eq. 7, HW:Eq. 6 when C1=0 and C2=1)
-    scalarField zEff(max((zDir() & pCf) - groundMin - d + z0, z0));
-
-    return sqr(Ustar(z0))/sqrt(Cmu_)*sqrt(C1_*log(zEff/z0) + C2_);
+    return
+        sqr(Ustar(z0))/sqrt(Cmu_)
+       *sqrt(C1_*log(((zDir() & pCf) - groundMin - d + z0)/z0) + C2_);
 }
 
 
@@ -251,9 +249,9 @@ tmp<scalarField> atmBoundaryLayer::epsilon(const vectorField& pCf) const
     const scalar groundMin = zDir() & ppMin_;
 
     // (YGCJ:Eq. 22; RH:Eq. 8, HW:Eq. 7 when C1=0 and C2=1)
-    scalarField zEff(max((zDir() & pCf) - groundMin - d + z0, z0));
-
-    return pow3(Ustar(z0))/(kappa_*zEff)*sqrt(C1_*log(zEff/z0) + C2_);
+    return
+        pow3(Ustar(z0))/(kappa_*((zDir() & pCf) - groundMin - d + z0))
+       *sqrt(C1_*log(((zDir() & pCf) - groundMin - d + z0)/z0) + C2_);
 }
 
 
@@ -265,9 +263,7 @@ tmp<scalarField> atmBoundaryLayer::omega(const vectorField& pCf) const
     const scalar groundMin = zDir() & ppMin_;
 
     // (YGJ:Eq. 13)
-    scalarField zEff(max((zDir() & pCf) - groundMin - d + z0, z0));
-
-    return Ustar(z0)/(kappa_*sqrt(Cmu_)*zEff);
+    return Ustar(z0)/(kappa_*sqrt(Cmu_)*((zDir() & pCf) - groundMin - d + z0));
 }
 
 

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

@@ -1,4 +1,4 @@
-  /*---------------------------------------------------------------------------*\
+/*---------------------------------------------------------------------------*\
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
@@ -36,18 +36,23 @@ License
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 template<class Type, class Stencil>
-void Foam::fv::LeastSquaresGrad<Type, Stencil>::calcGrad
-(
-    GeometricField
+Foam::tmp
+<
+    Foam::GeometricField
     <
-        typename outerProduct<vector, Type>::type,
-        fvPatchField,
-        volMesh
-    >& lsGrad,
-    const GeometricField<Type, fvPatchField, volMesh>& vtf
+        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();
 
@@ -57,7 +62,24 @@ void Foam::fv::LeastSquaresGrad<Type, Stencil>::calcGrad
         mesh
     );
 
-    lsGrad = dimensioned<GradType>(vtf.dimensions()/dimLength, Zero);
+    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();
     Field<GradType>& lsGradIf = lsGrad;
 
     const extendedCentredCellToCellStencil& stencil = lsv.stencil();
@@ -109,50 +131,6 @@ void 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,14 +132,6 @@ 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,21 +38,19 @@ License
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 template<class Type>
-void Foam::fv::fourthGrad<Type>::calcGrad
+Foam::tmp
+<
+    Foam::GeometricField
+    <
+        typename Foam::outerProduct<Foam::vector, Type>::type,
+        Foam::fvPatchField,
+        Foam::volMesh
+    >
+>
+Foam::fv::fourthGrad<Type>::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 word& name
 ) const
 {
     // The fourth-order gradient is calculated in two passes.  First,
@@ -61,11 +59,37 @@ void 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();
 
-    // Initialise to gradient
-    fGrad = secondfGrad;
+    // 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();
 
     const vectorField& C = mesh.C();
 
@@ -134,102 +158,9 @@ void 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,32 +99,6 @@ 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,67 +140,6 @@ 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,14 +145,6 @@ 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().topoChanging())
+    if (!this->mesh().cache(name) || this->mesh().changing())
     {
         // Delete any old occurrences to avoid double registration
         if (pgGrad && pgGrad->ownedByRegistry())
@@ -119,14 +119,17 @@ Foam::fv::gradScheme<Type>::grad
     }
     else
     {
-        if (pgGrad->upToDate(vsf) && this->mesh().upToDatePoints(*pgGrad))
+        if (pgGrad->upToDate(vsf))
         {
             solution::cachePrintMessage("Reusing", name, vsf);
         }
         else
         {
             solution::cachePrintMessage("Updating", name, vsf);
-            calcGrad(*pgGrad, vsf);
+            delete pgGrad;
+
+            pgGrad = calcGrad(vsf, name).ptr();
+            regIOobject::store(pgGrad);
         }
     }
 
@@ -177,23 +180,4 @@ 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,14 +142,6 @@ 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,25 +62,40 @@ void Foam::fv::cellLimitedGrad<Type, Limiter>::limitGradient
 
 
 template<class Type, class Limiter>
-void Foam::fv::cellLimitedGrad<Type, Limiter>::calcGrad
+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
+{
+    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,
-    const GeometricField<Type, fvPatchField, volMesh>& vsf
-) const
-{
-    const fvMesh& mesh = vsf.mesh();
-
-    basicGradScheme_().calcGrad(g, vsf);
-
-    if (k_ < SMALL)
-    {
-        return;
-    }
+    >& g = tGrad.ref();
 
     const labelUList& owner = mesh.owner();
     const labelUList& neighbour = mesh.neighbour();
@@ -212,49 +227,6 @@ void 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,14 +150,6 @@ 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,77 +2984,6 @@ 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,8 +53,6 @@ SourceFiles
 #include "lduPrimitiveMeshAssembly.H"
 #include "lduMesh.H"
 
-#include "fvMatrixExpression.H"
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -743,47 +741,6 @@ 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/fused/finiteVolume/fusedGaussGrad.C

@@ -226,94 +226,6 @@ 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,14 +145,6 @@ 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,87 +45,6 @@ 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,17 +65,6 @@ 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
         (
@@ -209,8 +198,6 @@ defineFvmLaplacianScalarGamma(sphericalTensor);
 defineFvmLaplacianScalarGamma(symmTensor);
 defineFvmLaplacianScalarGamma(tensor);
 
-#undef defineFvmLaplacianScalarGamma
-
 // Unweighted laplacian
 template<>
 tmp<GeometricField<scalar, fvPatchField, volMesh>>
@@ -240,6 +227,20 @@ 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,7 +28,6 @@ License
 
 #include "fusedGaussLaplacianScheme.H"
 #include "fvMesh.H"
-#include "fvMatrixExpression.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -40,83 +39,63 @@ template<>                                                                     \
 Foam::tmp<Foam::fvMatrix<Foam::Type>>                                          \
 Foam::fv::fusedGaussLaplacianScheme<Foam::Type, Foam::scalar>::                \
 fvmLaplacian                                                                   \
-(                                                                              \
-    const GeometricField<scalar, fvPatchField, volMesh>& gamma,                \
-    const GeometricField<Type, fvPatchField, volMesh>& vf                      \
-)                                                                              \
-{                                                                              \
-    DebugPout<< "fusedGaussLaplacianScheme::fvmLaplacian on " << vf.name()     \
-        << " with scalar gamma " << gamma.name() << " and ET" << endl;         \
-                                                                               \
-    const fvMesh& mesh = this->mesh();                                         \
-                                                                               \
-    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                                                   \
-    (                                                                          \
-        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::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;   \
+        << " with scalar gamma " << gamma.name() << 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                                                   \
+    GeometricField<scalar, fvsPatchField, surfaceMesh> gammaMagSf              \
     (                                                                          \
-        new fvMatrix<Type>                                                     \
-        (                                                                      \
-            vf,                                                                \
-            gamma.dimensions()*mesh.magSf().dimensions()*vf.dimensions()       \
-        )                                                                      \
+        gamma*mesh.magSf()                                                     \
+    );                                                                         \
+                                                                               \
+    tmp<fvMatrix<Type>> tfvm = fvmLaplacianUncorrected                         \
+    (                                                                          \
+        gammaMagSf,                                                            \
+        this->tsnGradScheme_().deltaCoeffs(vf),                                \
+        vf                                                                     \
     );                                                                         \
     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);              \
+        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();                                          \
+        }                                                                      \
     }                                                                          \
                                                                                \
     return tfvm;                                                               \
 }                                                                              \
                                                                                \
+                                                                               \
 template<>                                                                     \
 Foam::tmp<Foam::GeometricField<Foam::Type, Foam::fvPatchField, Foam::volMesh>> \
 Foam::fv::fusedGaussLaplacianScheme<Foam::Type, Foam::scalar>::                \
@@ -521,6 +500,38 @@ 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

@@ -120,7 +120,7 @@ dynamicContactAngleForce::dynamicContactAngleForce
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-dynamicContactAngleForce::~dynamicContactAngleForce()
+Foam::dynamicContactAngleForce::~dynamicContactAngleForce()
 {}  // distributionModel was forward declared
 
 

tutorials/incompressible/simpleFoam/pitzDaily_fused/Allrun

@@ -19,12 +19,4 @@ 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,15 +16,6 @@ FoamFile
 
 libs            (fusedFiniteVolume);
 
-/*
-DebugSwitches
-{
-    fusedGauss  1;
-    Gauss       1;
-    solution    1;
-}
-*/
-
 application     simpleFoam;
 
 startFrom       startTime;

tutorials/incompressible/simpleFoam/pitzDaily_fused/system/fvSolution

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

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

@@ -1,65 +0,0 @@
-/*--------------------------------*- 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;
-}
-
-
-// ************************************************************************* //