STYLE: ListExpression: copyright, indenting

Feature grad update

See merge request Development/openfoam!764

applications/test/expressionTemplates/Make/files

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

applications/test/expressionTemplates/Make/options

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

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

@@ -0,0 +1,126 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2025 M. Janssens
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    Test-expressionTemplates
+
+\*---------------------------------------------------------------------------*/
+
+#include "Time.H"
+#include "argList.H"
+#include "polyMesh.H"
+#include "pointMesh.H"
+#include "pointFields.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+
+    #include "createTime.H"
+    #include "createPolyMesh.H"
+
+    {
+        scalarField vals0({1.0, 2.0, 3.0});
+        const Expression::ListConstRefWrap<scalar> wvals0(vals0);
+
+        scalarField vals1;
+        Expression::ListRefWrap<scalar> wvals1(vals1.size(), vals1);
+        wvals1 = wvals0;
+        return 0;
+    }
+
+
+
+
+    const pointMesh& pMesh = pointMesh::New(mesh);
+
+    // Field, dimensionedScalar as List expression
+    {
+        scalarField vals({1.0, 2.0, 3.0});
+        dimensionedScalar d(dimless, 4.0);
+
+        scalarField result(d.expr(vals.size())*vals.expr());
+
+        Pout<< "result:" << result << endl;
+    }
+
+    Info<< "Reading field p\n" << endl;
+    pointScalarField p
+    (
+        IOobject
+        (
+            "pointDisplacement",
+            runTime.timeName(),
+            pMesh.thisDb(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        pMesh
+    );
+
+    pointScalarField result
+    (
+        IOobject
+        (
+            "result",
+            runTime.timeName(),
+            pMesh.thisDb(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            IOobject::NO_REGISTER
+        ),
+        pMesh,
+        dimensionedScalar("zero", Foam::sqr(p().dimensions()), 0)
+    );
+
+    // Construct value + dimensions with correct sizing
+    const auto oneField(dimensionedScalar(dimless, 1.0).expr(p));
+
+    // Make expression
+    auto expression = oneField * p.expr() + p.expr();
+
+    // Combine expressions
+    auto newExpression = sqr(expression);
+
+    // Assign values
+    result = newExpression;
+
+    DebugVar(result);
+
+
+    // Construct from expression
+    pointScalarField result2("result2", pMesh, sqrt(p.expr()));
+
+    DebugVar(result2);
+
+    return 0;
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,451 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    volPointInterpolationTest
+
+\*---------------------------------------------------------------------------*/
+
+#include "Time.H"
+#include "argList.H"
+#include "fvMesh.H"
+#include "ListExpression.H"
+#include "GeometricFieldExpression.H"
+#include "fvCFD.H"
+#include "fvMatrixExpression.H"
+#include <ratio>
+#include <chrono>
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+tmp<volScalarField> someFunction(const volScalarField& fld)
+{
+    return fld*1.0001;
+}
+
+
+template<class Type>
+void fusedGaussFvmLaplacian
+(
+    fvMatrix<Type>& fvm,
+    const surfaceInterpolationScheme<scalar>& interpGammaScheme,
+    const fv::snGradScheme<Type>& snGradScheme,
+    const GeometricField<scalar, fvPatchField, volMesh>& gamma,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    // Replacement for gaussLaplacianScheme::fvmLaplacian with scalar gamma
+    typedef GeometricField<Type, fvsPatchField, surfaceMesh> surfaceType;
+
+    const auto& mesh = vf.mesh();
+
+    // Expression for weights
+    const auto weights = interpGammaScheme.weights(gamma).expr();
+
+    // Expression for gamma_face * magSf
+    const auto gammaMagSf =
+        Expression::interpolate(gamma.expr(), weights, mesh)
+      * mesh.magSf().expr();
+
+    // Expression for deltaCoeffs
+    const auto deltaCoeffs = snGradScheme.deltaCoeffs(vf).expr();
+
+    // Construct matrix
+    Expression::fvmLaplacianUncorrected(fvm, gammaMagSf, deltaCoeffs);
+
+    if (snGradScheme.corrected())
+    {
+        // Wrap correction
+        const auto corr(snGradScheme.correction(vf).expr());
+        const auto V = mesh.V().expr();
+
+        if (mesh.fluxRequired(vf.name()))
+        {
+            fvm.faceFluxCorrectionPtr() = std::make_unique<surfaceType>
+            (
+                IOobject
+                (
+                    "faceFluxCorr",
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    IOobject::NO_REGISTER
+                ),
+                mesh,
+                gamma.dimensions()
+               *mesh.magSf().dimensions()
+               *corr.data().dimensions()
+            );
+            auto& faceFluxCorr = *fvm.faceFluxCorrectionPtr();
+            faceFluxCorr = gammaMagSf*corr;
+
+            fvm.source() =
+                fvm.source().expr()
+              - (
+                    V * fvc::div
+                    (
+                        faceFluxCorr
+                    )().primitiveField().expr()
+                );
+        }
+        else
+        {
+            // Temporary field
+            surfaceType faceFluxCorr
+            (
+                IOobject
+                (
+                    "faceFluxCorr",
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    IOobject::NO_REGISTER
+                ),
+                mesh,
+                gamma.dimensions()
+               *mesh.magSf().dimensions()
+               *corr.data().dimensions()
+            );
+            faceFluxCorr = gammaMagSf*corr;
+
+            fvm.source() =
+                fvm.source().expr()
+              - (
+                    V * fvc::div
+                    (
+                        faceFluxCorr
+                    )().primitiveField().expr()
+                );
+        }
+    }
+}
+
+
+using namespace std::chrono;
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh.thisDb(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+    
+    {
+        volScalarField p2("p2", p);
+        DebugVar(p2.boundaryFieldRef());
+    
+        for (auto& pf : p.boundaryFieldRef())
+        {
+            scalarField newVals(pf.size());
+            forAll(pf, i)
+            {
+                newVals[i] = scalar(i);
+            }
+            pf == newVals;
+        }
+
+        //std::vector<const scalarField*> ptrs;
+        // List<const scalarField*> ptrs;
+        // for (const auto& pp : p.boundaryField())
+        // {
+        //     ptrs.push_back(&pp);
+        // }
+        DebugVar(sizeof(long));
+        DebugVar(p.boundaryField());
+        Expression::ListsConstRefWrap<scalarField> expr(p.boundaryField());
+
+        const auto twoA = expr + expr;
+        Expression::ListsRefWrap<scalarField> expr2(p2.boundaryFieldRef());
+        Pout<< "**before assignment twoA:" << twoA.size()
+            << " expr2:" << expr2.size() << endl;
+        expr2 = twoA;
+        Pout<< "**after assignment twoA:" << twoA.size()
+            << " expr2:" << expr2.size() << endl;
+        DebugVar(p2.boundaryField());
+        // forAll(expr, i)
+        // {
+        //     Pout<< "i:" << i
+        //         //<< " expr:" << expr[i]
+        //         //<< " twoA:" << twoA[i]
+        //         << " expr2:" << expr2[i]
+        //         << endl;
+        // }
+        return 0;
+    }
+
+    {
+        //DebugVar(linearInterpolate(p));
+        //auto tweights = linear<scalar>(mesh).weights(p);
+        //DebugVar(tweights);
+
+        surfaceScalarField result
+        (
+            IOobject
+            (
+                "result",
+                runTime.timeName(),
+                mesh.thisDb(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                IOobject::NO_REGISTER
+            ),
+            mesh,
+            dimensionedScalar(p.dimensions(), 0)
+        );
+        //result = Expression::interpolate
+        //(
+        //    p.expr(),
+        //    tweights().expr(),
+        //    mesh
+        //);
+
+        result = Expression::linearInterpolate(p.expr(), mesh);
+
+
+        DebugVar(result);
+
+        return 0;
+    }
+
+
+
+    volScalarField p2
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh.thisDb(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE,
+            IOobject::NO_REGISTER
+        ),
+        mesh
+    );
+
+    // Expresions of volFields
+    {
+        volScalarField result
+        (
+            IOobject
+            (
+                "result",
+                runTime.timeName(),
+                mesh.thisDb(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                IOobject::NO_REGISTER
+            ),
+            mesh,
+            dimensionedScalar(p.dimensions(), 0)
+        );
+
+        {
+            Pout<< "No expression templates:" << endl;
+            const high_resolution_clock::time_point t1 =
+                high_resolution_clock::now();
+            result = p + p2;
+            const high_resolution_clock::time_point t2 =
+                high_resolution_clock::now();
+            const duration<double> time_span = t2 - t1;
+            Pout<< "Operation time:" << time_span.count() << endl;
+        }
+        {
+            Pout<< "With expression templates:" << endl;
+            const high_resolution_clock::time_point t1 =
+                high_resolution_clock::now();
+            result = p.expr() + p2.expr();
+            const high_resolution_clock::time_point t2 =
+                high_resolution_clock::now();
+            const duration<double> time_span = t2 - t1;
+            Pout<< "Operation time:" << time_span.count() << endl;
+        }
+
+//        const auto oldDimensions = p.dimensions();
+//        p.dimensions().reset(dimless);
+//        p2.dimensions().reset(dimless);
+//        result.dimensions().reset(dimless);
+//        {
+//
+//            Pout<< "Complex expression : No expression templates:" << endl;
+//            const high_resolution_clock::time_point t1 =
+//                high_resolution_clock::now();
+//            result = cos(p + 0.5*sqrt(p2-sin(p)));
+//            const high_resolution_clock::time_point t2 =
+//                high_resolution_clock::now();
+//            const duration<double> time_span = t2 - t1;
+//            Pout<< "Operation time:" << time_span.count() << endl;
+//        }
+//        {
+//            Pout<< "Complex expression : With expression templates:" << endl;
+//            const high_resolution_clock::time_point t1 =
+//                high_resolution_clock::now();
+//            const auto zeroDotFive
+//            (
+//                dimensionedScalar(dimless, 0.5).expr(p)
+//            );
+//            result = cos(p.expr() + zeroDotFive*sqrt(p2.expr()-sin(p.expr())));
+//            const high_resolution_clock::time_point t2 =
+//                high_resolution_clock::now();
+//            const duration<double> time_span = t2 - t1;
+//            Pout<< "Operation time:" << time_span.count() << endl;
+//        }
+//        p.dimensions().reset(oldDimensions);
+//        p2.dimensions().reset(oldDimensions);
+//        result.dimensions().reset(oldDimensions);
+
+        return 0;
+//        auto expression = someFunction(p).expr() + someFunction(p).expr();
+//        result = expression;
+//        DebugVar(result);
+    }
+    {
+        volScalarField result
+        (
+            IOobject
+            (
+                "result",
+                runTime.timeName(),
+                mesh.thisDb(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                IOobject::NO_REGISTER
+            ),
+            mesh,
+            dimensionedScalar(p.dimensions(), 0)
+        );
+        auto expression = someFunction(p).expr() + someFunction(p).expr();
+        result = expression;
+        DebugVar(result);
+    }
+
+    // Expresions of volFields
+    {
+        volScalarField result
+        (
+            "result",
+            mesh,
+            sqr(p.expr() + p.expr())
+        );
+        DebugVar(result);
+    }
+    {
+        // Fill p with some values
+        forAll(p, celli)
+        {
+            p[celli] = celli;
+        }
+        p.correctBoundaryConditions();
+
+        // Interpolate to surface field
+        surfaceScalarField result
+        (
+            "result",
+            mesh,
+            Expression::interpolate    //<volExpr, surfaceExpr>
+            (
+                p.expr(),
+                mesh.surfaceInterpolation::weights().expr(),
+                mesh
+            )
+        );
+        DebugVar(result);
+    }
+    {
+        // For testing as a replacement of laplacian weights
+        const volScalarField gamma
+        (
+            IOobject
+            (
+                "gamma",
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                IOobject::NO_REGISTER
+            ),
+            mesh,
+            dimensionedScalar(dimless, 1.0)
+        );
+
+        fvMatrix<scalar> fvm
+        (
+            p,
+            gamma.dimensions()*mesh.magSf().dimensions()*p.dimensions()
+        );
+
+        const linear<scalar> interpGammaScheme(mesh);
+        const fv::correctedSnGrad<scalar> snGradScheme(mesh);
+
+        fusedGaussFvmLaplacian
+        (
+            fvm,
+            interpGammaScheme,
+            snGradScheme,
+            gamma,
+            p
+        );
+
+        DebugVar(fvm.source());
+    }
+
+
+
+    // Expressions of fvMatrix
+    {
+        tmp<fvMatrix<scalar>> tm0(fvm::laplacian(p));
+        const fvMatrix<scalar>& m0 = tm0();
+        DebugVar(m0.dimensions());
+
+        tmp<fvMatrix<scalar>> tm1(fvm::laplacian(p));
+        const fvMatrix<scalar>& m1 = tm1();
+        DebugVar(m1.dimensions());
+
+        fvMatrix<scalar> m2(p, m0.expr() + m1.expr());
+        DebugVar(m2.dimensions());
+    }
+
+    return 0;
+}
+
+
+// ************************************************************************* //

src/Allwmake

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

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

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

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

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

src/OpenFOAM/dimensionedTypes/dimensionedType/dimensionedType.H

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

src/OpenFOAM/expressionTemplates/GenericExpression.H

@@ -0,0 +1,252 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2025 M. Janssens
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Namespace
+    Foam::GenericExpression
+
+Description
+    Expression templates.
+
+SourceFiles
+    GenericExpression.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_GenericExpression_H
+#define Foam_GenericExpression_H
+
+#include <cassert>
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace Expression
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class GenericExpression Declarations
+\*---------------------------------------------------------------------------*/
+
+template<typename E>
+class GenericExpression
+{
+public:
+    static constexpr bool is_leaf = false;
+
+    auto evaluate() const
+    {
+        // Delegation to the actual expression type. This avoids dynamic
+        // polymorphism (a.k.a. virtual functions in C++)
+        return static_cast<E const&>(*this).evaluate();
+    }
+};
+
+
+// Macros
+// ~~~~~~
+
+#define EXPRESSION_FUNCTION1(BaseClass, Func, BaseFunc, OpFunc)                \
+template<typename E1>                                                          \
+class OpFunc                                                                   \
+:                                                                              \
+    public BaseClass<OpFunc<E1>>                                               \
+{                                                                              \
+    typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_;      \
+                                                                               \
+public:                                                                        \
+    static constexpr bool is_leaf = false;                                     \
+                                                                               \
+    OpFunc(const E1& u)                                                        \
+    :                                                                          \
+        u_(u)                                                                  \
+    {}                                                                         \
+                                                                               \
+    auto evaluate() const                                                      \
+    {                                                                          \
+        return BaseFunc(u_.evaluate());                                        \
+    }                                                                          \
+};                                                                             \
+template<typename E1>                                                          \
+OpFunc<E1> Func                                                                \
+(                                                                              \
+    const BaseClass<E1>& u                                                     \
+)                                                                              \
+{                                                                              \
+    return OpFunc<E1>(static_cast<const E1&>(u));                              \
+}
+
+
+#define EXPRESSION_FUNCTION2(BaseClass, Func, BaseFunc, OpFunc)                \
+template<typename E1, typename E2>                                             \
+class OpFunc                                                                   \
+:                                                                              \
+    public BaseClass<OpFunc<E1, E2>>                                           \
+{                                                                              \
+    typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_;      \
+    typename std::conditional<E2::is_leaf, const E2&, const E2>::type v_;      \
+                                                                               \
+public:                                                                        \
+    static constexpr bool is_leaf = false;                                     \
+                                                                               \
+    OpFunc(const E1& u, const E2& v)                                           \
+    :                                                                          \
+     u_(u), v_(v)                                                              \
+    {}                                                                         \
+    auto evaluate() const                                                      \
+    {                                                                          \
+        return BaseFunc(u_.evaluate(), v_.evaluate());                         \
+    }                                                                          \
+};                                                                             \
+template<typename E1, typename E2>                                             \
+OpFunc<E1, E2> Func                                                            \
+(                                                                              \
+    const BaseClass<E1>& u,                                                    \
+    const BaseClass<E2>& v                                                     \
+)                                                                              \
+{                                                                              \
+    return OpFunc<E1, E2>                                                      \
+    (                                                                          \
+        static_cast<const E1&>(u),                                             \
+        static_cast<const E2&>(v)                                              \
+    );                                                                         \
+}
+
+
+#define EXPRESSION_OPERATOR(BaseClass, Op, BaseOp, OpFunc)                     \
+template<typename E1, typename E2>                                             \
+class OpFunc                                                                   \
+:                                                                              \
+    public BaseClass                                                           \
+    <                                                                          \
+        OpFunc<E1, E2>                                                         \
+    >                                                                          \
+{                                                                              \
+    /* cref if leaf, copy otherwise */                                         \
+    typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_;      \
+    typename std::conditional<E2::is_leaf, const E2&, const E2>::type v_;      \
+                                                                               \
+public:                                                                        \
+    static constexpr bool is_leaf = false;                                     \
+                                                                               \
+    OpFunc(const E1& u, const E2& v)                                           \
+    :                                                                          \
+     u_(u), v_(v)                                                              \
+    {}                                                                         \
+    auto evaluate() const                                                      \
+    {                                                                          \
+        return u_.evaluate() BaseOp v_.evaluate();                             \
+    }                                                                          \
+};                                                                             \
+template<typename E1, typename E2>                                             \
+OpFunc<E1, E2>                                                                 \
+operator Op                                                                    \
+(                                                                              \
+    BaseClass<E1> const& u,                                                    \
+    BaseClass<E2> const& v                                                     \
+)                                                                              \
+{                                                                              \
+    return OpFunc<E1, E2>                                                      \
+    (                                                                          \
+        static_cast<const E1&>(u),                                             \
+        static_cast<const E2&>(v)                                              \
+    );                                                                         \
+}
+
+
+// Do '-' separately until we work out macro expansion...
+template<typename E1>
+class g_negate
+:
+    public GenericExpression
+    <
+        g_negate<E1>
+    >
+{
+    typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_;
+
+public:
+    static constexpr bool is_leaf = false;
+
+    g_negate(const E1& u)
+    :
+        u_(u)
+    {}
+
+    auto evaluate() const
+    {
+        return -u_.evaluate();
+    }
+};
+template<typename E1>
+g_negate<E1> operator-
+(
+    const GenericExpression<E1>& u
+)
+{
+    return g_negate<E1>(static_cast<const E1&>(u));
+}
+
+
+EXPRESSION_FUNCTION1(GenericExpression, sqr, Foam::sqr, g_sqr)
+EXPRESSION_FUNCTION1(GenericExpression, sqrt, Foam::sqrt, g_sqrt)
+EXPRESSION_FUNCTION1(GenericExpression, magSqr, Foam::magSqr, g_magSqr)
+EXPRESSION_FUNCTION1(GenericExpression, symm, Foam::symm, g_symm)
+EXPRESSION_FUNCTION1(GenericExpression, pow2, Foam::pow2, g_pow2)
+EXPRESSION_FUNCTION1(GenericExpression, pow3, Foam::pow3, g_pow3)
+EXPRESSION_FUNCTION1(GenericExpression, pow4, Foam::pow4, g_pow4)
+EXPRESSION_FUNCTION1(GenericExpression, operator~, Foam::operator~, g_compl)
+
+//TBD. Parse problem
+//EXPRESSION_FUNCTION1(GenericExpression, operator!, Foam::operator!, g_not)
+
+#undef EXPRESSION_FUNCTION1
+
+EXPRESSION_FUNCTION2(GenericExpression, operator+, Foam::operator+, g_add);
+EXPRESSION_FUNCTION2(GenericExpression, operator-, Foam::operator-, g_subtract);
+EXPRESSION_FUNCTION2(GenericExpression, operator*, Foam::operator*, g_multiply);
+EXPRESSION_FUNCTION2(GenericExpression, operator/, Foam::operator/, g_divide);
+
+#undef EXPRESSION_FUNCTION2
+
+EXPRESSION_OPERATOR(GenericExpression, ||, ||, g_or)
+EXPRESSION_OPERATOR(GenericExpression, &&, &&, g_and)
+EXPRESSION_OPERATOR(GenericExpression, &, &, g_bitand)
+EXPRESSION_OPERATOR(GenericExpression, |, |, g_bitor)
+EXPRESSION_OPERATOR(GenericExpression, ^, ^, g_xor)
+
+#undef EXPRESSION_OPERATOR
+
+} // End namespace Expression
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/expressionTemplates/boolExpression.H

@@ -0,0 +1,145 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2025 M. Janssens
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Namespace
+    Foam::boolExpression
+
+Description
+    Expression templates for bool
+
+SourceFiles
+    boolExpression.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_boolExpression_H
+#define Foam_boolExpression_H
+
+#include "GenericExpression.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace Expression
+{
+
+/*---------------------------------------------------------------------------*\
+                       Namespace Expression Declarations
+\*---------------------------------------------------------------------------*/
+
+// Wrap of non-const bool
+// ~~~~~~~~~~~~~~~~~~~~~~
+
+class boolWrap
+:
+    public GenericExpression<boolWrap>
+{
+    bool elems_;
+
+public:
+    static constexpr bool is_leaf = false;  //true;
+
+    //- Construct from components
+    boolWrap(const bool elems)
+    :
+        elems_(elems)
+    {}
+
+    // Construct from GenericExpression, forcing its evaluation.
+    template<typename E>
+    boolWrap
+    (
+        const GenericExpression<E>& expr
+    )
+    :
+        elems_(expr.evaluate())
+    {}
+
+    auto evaluate() const
+    {
+        return elems_;
+    }
+
+    template<typename E>
+    auto evaluate
+    (
+        const GenericExpression<E>& expr
+    )
+    {
+        elems_ = expr.evaluate();
+        return elems_;
+    }
+
+    //- Assignment
+    template<typename E>
+    void operator=
+    (
+        const GenericExpression<E>& expr
+    )
+    {
+        elems_ = expr.evaluate();
+    }
+};
+
+
+// Wrap of const bool
+// ~~~~~~~~~~~~~~~~~~
+
+class boolConstWrap
+:
+    public GenericExpression<boolConstWrap>
+{
+    const bool elems_;
+
+public:
+
+    // ! Store as copy (since holds reference)
+    static constexpr bool is_leaf = false;
+
+    // construct from components
+    boolConstWrap(const bool elems)
+    :
+        elems_(elems)
+    {}
+
+    auto evaluate() const
+    {
+        return elems_;
+    }
+};
+
+
+} // End namespace Expression
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/expressionTemplates/dimensionSetExpression.H

@@ -0,0 +1,148 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2025 M. Janssens
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Namespace
+    Foam::dimensionSetExpression
+
+Description
+    Expression templates for dimensionSet
+
+SourceFiles
+    dimensionSetExpression.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_dimensionSetExpression_H
+#define Foam_dimensionSetExpression_H
+
+#include "GenericExpression.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace Expression
+{
+
+/*---------------------------------------------------------------------------*\
+                       Namespace Expression Declarations
+\*---------------------------------------------------------------------------*/
+
+// Wrap of non-const reference to dimensionSet
+// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+class dimensionSetRefWrap
+:
+    public GenericExpression<dimensionSetRefWrap>
+{
+    dimensionSet& elems_;
+
+public:
+    static constexpr bool is_leaf = false;  //true;
+
+    //- Construct from components
+    dimensionSetRefWrap(dimensionSet& elems)
+    :
+        elems_(elems)
+    {}
+
+    // Construct from GenericExpression, forcing its evaluation.
+    template<typename E>
+    dimensionSetRefWrap
+    (
+        dimensionSet& elems,
+        const GenericExpression<E>& expr
+    )
+    :
+        elems_(elems)
+    {
+        elems_ = expr.evaluate();
+    }
+
+    auto evaluate() const
+    {
+        return elems_;
+    }
+
+    template<typename E>
+    auto& evaluate
+    (
+        const GenericExpression<E>& expr
+    )
+    {
+        elems_ = expr.evaluate();
+        return elems_;
+    }
+
+    //- Assignment
+    template<typename E>
+    void operator=
+    (
+        const GenericExpression<E>& expr
+    )
+    {
+        elems_ = expr.evaluate();
+    }
+};
+
+
+// Wrap of const reference to dimensionSet
+// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+class dimensionSetConstRefWrap
+:
+    public GenericExpression<dimensionSetConstRefWrap>
+{
+    const dimensionSet& elems_;
+
+public:
+
+    // ! Store as copy (since holds reference)
+    static constexpr bool is_leaf = false;
+
+    // construct from components
+    dimensionSetConstRefWrap(const dimensionSet& elems)
+    :
+        elems_(elems)
+    {}
+
+    auto evaluate() const
+    {
+        return elems_;
+    }
+};
+
+
+} // End namespace Expression
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

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

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

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

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

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

src/OpenFOAM/memory/tmp/tmp.H

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

src/OpenFOAM/memory/tmp/tmpI.H

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C

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

src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H

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

src/fused/finiteVolume/fusedGaussGrad.C

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

src/fused/finiteVolume/fusedGaussGrad.H

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

src/fused/finiteVolume/fusedGaussLaplacianScheme.C

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

src/fused/finiteVolume/fusedGaussLaplacianScheme.H

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

src/fused/finiteVolume/fusedGaussLaplacianSchemes.C

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

src/regionFaModels/Make/files

@@ -25,20 +25,11 @@ $(kinematic)/injectionModel/injectionModelList/injectionModelList.C
 $(kinematic)/injectionModel/injectionModel/injectionModel.C
 $(kinematic)/injectionModel/injectionModel/injectionModelNew.C
 
-/* Film separation models */
-
-filmSeparation=$(kinematic)/injectionModel/filmSeparation
-filmSeparationModels=$(filmSeparation)/filmSeparationModels
-$(filmSeparation)/filmSeparation.C
-$(filmSeparationModels)/filmSeparationModel/filmSeparationModel.C
-$(filmSeparationModels)/filmSeparationModel/filmSeparationModelNew.C
-$(filmSeparationModels)/OwenRyleyModel/OwenRyleyModel.C
-$(filmSeparationModels)/FriedrichModel/FriedrichModel.C
-cornerDetectionModels=$(filmSeparation)/cornerDetectionModels
-$(cornerDetectionModels)/cornerDetectionModel/cornerDetectionModel.C
-$(cornerDetectionModels)/cornerDetectionModel/cornerDetectionModelNew.C
-$(cornerDetectionModels)/fluxBased/fluxBased.C
-$(cornerDetectionModels)/geometryBased/geometryBased.C
+$(kinematic)/injectionModel/filmSeparation/filmSeparation.C
+$(kinematic)/injectionModel/filmSeparation/filmSeparationModels/filmSeparationModel/filmSeparationModel.C
+$(kinematic)/injectionModel/filmSeparation/filmSeparationModels/filmSeparationModel/filmSeparationModelNew.C
+$(kinematic)/injectionModel/filmSeparation/filmSeparationModels/OwenRyleyModel/OwenRyleyModel.C
+$(kinematic)/injectionModel/filmSeparation/filmSeparationModels/FriedrichModel/FriedrichModel.C
 
 $(kinematic)/injectionModel/BrunDrippingInjection/BrunDrippingInjection.C
 

src/regionFaModels/Make/options

@@ -11,8 +11,7 @@ EXE_INC = \
     -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
-    -I$(LIB_SRC)/transportModels \
-    -I$(LIB_SRC)/fileFormats/lnInclude
+    -I$(LIB_SRC)/transportModels
 
 LIB_LIBS = \
     -lfiniteVolume \
@@ -25,5 +24,4 @@ LIB_LIBS = \
     -lthermophysicalProperties \
     -lspecie \
     -lfaOptions \
-    -ldistributionModels \
-    -lfileFormats
+    -ldistributionModels

src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/cornerDetectionModels/cornerDetectionModel/cornerDetectionModel.C

@@ -1,103 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2025 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "cornerDetectionModel.H"
-#include "faMesh.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTypeNameAndDebug(cornerDetectionModel, 0);
-    defineRunTimeSelectionTable(cornerDetectionModel, dictionary);
-}
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-Foam::scalar Foam::cornerDetectionModel::dihedralAngle
-(
-    const vector& n0,
-    const vector& n1
-) const
-{
-    if (mag(n0) <= VSMALL || mag(n1) <= VSMALL)
-    {
-#ifdef FULL_DEBUG
-        WarningInFunction
-            << "Degenerate face normal magnitude (|n| ~ 0). "
-            << "Returning 0 for dihedral angle." << nl;
-#endif
-        return 0;
-    }
-
-    const scalar a = mag(n1 - n0);
-    const scalar b = mag(n1 + n0);
-
-    // The dihedral angle is calculated as 2*atan2(|n1 - n0|, |n1 + n0|),
-    // which gives the angle between the two normals n0 and n1.
-    scalar phi = scalar(2)*std::atan2(a, b);
-
-    // Clamp to [0, pi]
-    phi = max(0, min(constant::mathematical::pi, phi));
-
-    if (!std::isfinite(phi))
-    {
-#ifdef FULL_DEBUG
-        WarningInFunction
-            << "Non-finite dihedral angle computed. Returning 0." << nl;
-#endif
-        return 0;
-    }
-
-    return phi;
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::cornerDetectionModel::cornerDetectionModel
-(
-    const faMesh& mesh,
-    const regionModels::areaSurfaceFilmModels::liquidFilmBase& film,
-    const dictionary& dict
-)
-:
-    mesh_(mesh),
-    film_(film),
-    dict_(dict)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::cornerDetectionModel::~cornerDetectionModel()
-{}  // faMesh was forward declared
-
-
-// ************************************************************************* //

src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/cornerDetectionModels/cornerDetectionModel/cornerDetectionModel.H

@@ -1,211 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2025 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Namespace
-    Foam::cornerDetectionModels
-
-Description
-    A namespace for various \c cornerDetection model implementations.
-
-Class
-    Foam::cornerDetectionModel
-
-Description
-    A base class for \c cornerDetection models.
-
-SourceFiles
-    cornerDetectionModel.C
-    cornerDetectionModelNew.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef Foam_cornerDetectionModel_H
-#define Foam_cornerDetectionModel_H
-
-#include "liquidFilmBase.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward Declarations
-class faMesh;
-class dictionary;
-
-/*---------------------------------------------------------------------------*\
-                    Class cornerDetectionModel Declaration
-\*---------------------------------------------------------------------------*/
-
-class cornerDetectionModel
-{
-    // Private Data
-
-        //- Const reference to the finite-area mesh
-        const faMesh& mesh_;
-
-        //- Const reference to the film
-        const regionModels::areaSurfaceFilmModels::liquidFilmBase& film_;
-
-        //- Const reference to the dictionary
-        const dictionary& dict_;
-
-        //- Bitset for corner faces, true for corner faces
-        bitSet cornerFaces_;
-
-        //- List of corner angles [rad]
-        scalarList cornerAngles_;
-
-
-protected:
-
-    // Protected Member Functions
-
-        //- Return the dihedral angle between two neighbour faces
-        //  Robust 2*atan form (handles parallel normals better than acos)
-        //  \param n0 First normal vector
-        //  \param n1 Second normal vector
-        //  \return The dihedral angle [rad]
-        scalar dihedralAngle(const vector& n0, const vector& n1) const;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("cornerDetectionModel");
-
-
-    // Declare runtime constructor selection table
-
-        declareRunTimeSelectionTable
-        (
-            autoPtr,
-            cornerDetectionModel,
-            dictionary,
-            (
-                const faMesh& mesh,
-                const regionModels::areaSurfaceFilmModels::liquidFilmBase& film,
-                const dictionary& dict
-            ),
-            (mesh, film, dict)
-        );
-
-
-    // Selectors
-
-        //- Return a reference to the selected cornerDetection model
-        static autoPtr<cornerDetectionModel> New
-        (
-            const faMesh& mesh,
-            const regionModels::areaSurfaceFilmModels::liquidFilmBase& film,
-            const dictionary& dict
-        );
-
-
-    // Constructors
-
-        //- Construct from components
-        cornerDetectionModel
-        (
-            const faMesh& mesh,
-            const regionModels::areaSurfaceFilmModels::liquidFilmBase& film,
-            const dictionary& dict
-        );
-
-
-    // Generated Methods
-
-        //- No copy construct
-        cornerDetectionModel(const cornerDetectionModel&) = delete;
-
-        //- No copy assignment
-        void operator=(const cornerDetectionModel&) = delete;
-
-
-    //- Destructor
-    virtual ~cornerDetectionModel();
-
-
-    // Member Functions
-
-    // Getters
-
-        //- Return const reference to the finite-area mesh
-        const faMesh& mesh() const noexcept { return mesh_; }
-
-        //- Return const reference to the film
-        const regionModels::areaSurfaceFilmModels::liquidFilmBase&
-            film() const noexcept { return film_; }
-
-        //- Return const reference to the dictionary
-        const dictionary& dict() const noexcept { return dict_; }
-
-        //- Return const reference to the corner faces bitset
-        const bitSet& getCornerFaces() const noexcept { return cornerFaces_; }
-
-        //- Return const reference to the corner angles list [rad]
-        const scalarList& getCornerAngles() const noexcept
-        {
-            return cornerAngles_;
-        }
-
-
-    // Setters
-
-        //- Set the corner faces bitset
-        void setCornerFaces(const bitSet& cornerFaces)
-        {
-            cornerFaces_ = cornerFaces;
-        }
-
-        //- Set the corner angles list [rad]
-        void setCornerAngles(const scalarList& cornerAngles)
-        {
-            cornerAngles_ = cornerAngles;
-        }
-
-
-    // Evaluation
-
-        //- Detect and store the corner faces
-        virtual bool detectCorners() = 0;
-
-
-    // I-O
-
-        //- Read the model dictionary
-        virtual bool read(const dictionary&) { return true; }
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/cornerDetectionModels/cornerDetectionModel/cornerDetectionModelNew.C

@@ -1,65 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2025 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "cornerDetectionModel.H"
-#include "faMesh.H"
-#include "liquidFilmBase.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-Foam::autoPtr<Foam::cornerDetectionModel> Foam::cornerDetectionModel::New
-(
-    const faMesh& mesh,
-    const regionModels::areaSurfaceFilmModels::liquidFilmBase& film,
-    const dictionary& dict
-)
-{
-    const word modelType
-    (
-        dict.getOrDefault<word>("cornerDetectionModel", "fluxBased")
-    );
-
-    Info<< "    Selecting corner-detection model: " << modelType << nl << endl;
-
-    auto* ctorPtr = dictionaryConstructorTable(modelType);
-
-    if (!ctorPtr)
-    {
-        FatalIOErrorInLookup
-        (
-            dict,
-            "cornerDetectionModel",
-            modelType,
-            *dictionaryConstructorTablePtr_
-        ) << exit(FatalIOError);
-    }
-
-    return autoPtr<cornerDetectionModel>(ctorPtr(mesh, film, dict));
-}
-
-
-// ************************************************************************* //

src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/cornerDetectionModels/fluxBased/fluxBased.C

@@ -1,391 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2025 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fluxBased.H"
-#include "OBJstream.H"
-#include "processorFaPatch.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace cornerDetectionModels
-{
-    defineTypeNameAndDebug(fluxBased, 0);
-    addToRunTimeSelectionTable(cornerDetectionModel, fluxBased, dictionary);
-}
-}
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-Foam::bitSet Foam::cornerDetectionModels::fluxBased::identifyCornerEdges() const
-{
-    // Return true if face normals converge, i.e. sharp edge
-    // Face-normal vectors diverge: no separation, converge: separation (maybe)
-    const auto isCornerEdgeSharp =
-    [](
-        const vector& fcO,  // face-centre owner
-        const vector& fcN,  // face-centre neigh
-        const vector& fnO,  // face-normal owner
-        const vector& fnN   // face-normal neigh
-    ) noexcept -> bool
-    {
-        // Threshold for sharpness detection
-        constexpr scalar sharpEdgeThreshold = -1e-8;
-
-        // Relative centre and normal of the two faces sharing the edge
-        const vector relativePosition(fcN - fcO);
-        const vector relativeNormal(fnN - fnO);
-
-        // Sharp if normals converge along the centre-to-centre direction
-        return ((relativeNormal & relativePosition) < sharpEdgeThreshold);
-    };
-
-
-    // Cache the operand references
-    const areaVectorField& fc = mesh().areaCentres();
-    const areaVectorField& fn = mesh().faceAreaNormals();
-    const labelUList& own = mesh().edgeOwner();
-    const labelUList& nei = mesh().edgeNeighbour();
-
-
-    // Allocate the resource for the return object
-    bitSet cornerEdges(mesh().nEdges(), false);
-
-    // Internal edges (owner <-> neighbour)
-    forAll(nei, edgei)
-    {
-        const label faceO = own[edgei];
-        const label faceN = nei[edgei];
-
-        cornerEdges[edgei] = isCornerEdgeSharp
-        (
-            fc[faceO],
-            fc[faceN],
-            fn[faceO],
-            fn[faceN]
-        );
-    }
-
-
-    // Skip the rest of the routine if the simulation is a serial run
-    if (!Pstream::parRun()) return cornerEdges;
-
-    // Check if processor face-normal vectors diverge (no separation)
-    // or converge (separation may occur)
-    const faBoundaryMesh& patches = mesh().boundary();
-
-    for (const faPatch& fap : patches)
-    {
-        if (!isA<processorFaPatch>(fap)) continue;
-
-        const label patchi = fap.index();
-        const auto& edgeFaces = fap.edgeFaces();
-        const label internalEdge0 = fap.start();
-
-        const auto& fcp = fc.boundaryField()[patchi];
-        const auto& fnp = fn.boundaryField()[patchi];
-
-        // Processor edges (owner <-| none)
-        forAll(fnp, bndEdgei)
-        {
-            const label faceO = edgeFaces[bndEdgei];
-            const label meshEdgei = internalEdge0 + bndEdgei;
-
-            cornerEdges[meshEdgei] = isCornerEdgeSharp
-            (
-                fc[faceO],
-                fcp[bndEdgei],
-                fn[faceO],
-                fnp[bndEdgei]
-            );
-        }
-    }
-
-    return cornerEdges;
-}
-
-
-Foam::bitSet Foam::cornerDetectionModels::fluxBased::identifyCornerFaces
-(
-    const bitSet& cornerEdges
-) const
-{
-    // Marks the separating face based on edge flux sign
-    const auto markSeparation =
-    [](
-        bitSet& cornerFaces,
-        const scalar phiEdge,
-        const label faceO,
-        const label faceN = -1  /* = -1 for processor edges */
-    ) noexcept -> void
-    {
-        constexpr scalar tol = 1e-8;
-
-        // Assuming no sources/sinks at the edge
-        if (phiEdge > tol)  // From owner to neighbour
-        {
-            cornerFaces[faceO] = true;
-        }
-        else if ((phiEdge < -tol) && (faceN != -1))  // From nei to own
-        {
-            cornerFaces[faceN] = true;
-        }
-    };
-
-
-    // Cache the operand references
-    const edgeScalarField& phis = film().phi2s();
-    const labelUList& own = mesh().edgeOwner();
-    const labelUList& nei = mesh().edgeNeighbour();
-
-    // Allocate the resource for the return object
-    bitSet cornerFaces(mesh().faces().size(), false);
-
-    // Internal faces (owner <-> neighbour)
-    forAll(nei, edgei)
-    {
-        if (!cornerEdges[edgei]) continue;
-
-        markSeparation
-        (
-            cornerFaces,
-            phis[edgei],
-            own[edgei],  // faceO
-            nei[edgei]   // faceN
-        );
-    }
-
-
-    // Skip the rest of the routine if the simulation is a serial run
-    if (!Pstream::parRun()) return cornerFaces;
-
-    const faBoundaryMesh& patches = mesh().boundary();
-
-    for (const faPatch& fap : patches)
-    {
-        if (!isA<processorFaPatch>(fap)) continue;
-
-        const label patchi = fap.index();
-        const auto& edgeFaces = fap.edgeFaces();
-        const label internalEdge0 = fap.start();
-
-        const auto& phisp = phis.boundaryField()[patchi];
-
-        // Processor faces (owner <-| none)
-        forAll(phisp, bndEdgei)
-        {
-            const label faceO = edgeFaces[bndEdgei];
-            const label meshEdgei = internalEdge0 + bndEdgei;
-
-            if (!cornerEdges[meshEdgei]) continue;
-
-            markSeparation
-            (
-                cornerFaces,
-                phisp[bndEdgei],
-                faceO
-                /*faceN = -1*/
-            );
-        }
-    }
-
-    return cornerFaces;
-}
-
-
-Foam::scalarList Foam::cornerDetectionModels::fluxBased::calcCornerAngles
-(
-    const bitSet& faces,
-    const bitSet& edges
-) const
-{
-    // Cache the operand references
-    const areaVectorField& fn = mesh().faceAreaNormals();
-    const labelUList& own = mesh().edgeOwner();
-    const labelUList& nei = mesh().edgeNeighbour();
-
-    scalarList cornerFaceAngles(mesh().faces().size(), Zero);
-
-    // Internal edges (owner <-> neighbour)
-    forAll(nei, edgei)
-    {
-        if (!edges[edgei]) continue;
-
-        const label faceO = own[edgei];
-        const label faceN = nei[edgei];
-
-        // If neither adjacent face is flagged as a corner, skip the atan2 work
-        if (!faces[faceO] && !faces[faceN]) continue;
-
-        const scalar ang = this->dihedralAngle(fn[faceO], fn[faceN]);
-
-        if (faces[faceO]) cornerFaceAngles[faceO] = ang;
-        if (faces[faceN]) cornerFaceAngles[faceN] = ang;
-    }
-
-
-    // Skip the rest of the routine if the simulation is a serial run
-    if (!Pstream::parRun()) return cornerFaceAngles;
-
-    const faBoundaryMesh& patches = mesh().boundary();
-
-    for (const faPatch& fap : patches)
-    {
-        if (!isA<processorFaPatch>(fap)) continue;
-        const label patchi = fap.index();
-        const auto& edgeFaces = fap.edgeFaces();
-        const label internalEdge0 = fap.start();
-
-        const auto& fnp = fn.boundaryField()[patchi];
-
-        // Processor edges (owner <-| none)
-        forAll(fnp, bndEdgei)
-        {
-            const label faceO = edgeFaces[bndEdgei];
-            const label meshEdgei = internalEdge0 + bndEdgei;
-
-                // Only if the mesh edge and owner face are both corners
-                if (!edges[meshEdgei] || !faces[faceO]) continue;
-
-                cornerFaceAngles[faceO] =
-                    this->dihedralAngle(fn[faceO], fnp[bndEdgei]);
-        }
-    }
-
-    return cornerFaceAngles;
-}
-
-
-void Foam::cornerDetectionModels::fluxBased::writeEdgesAndFaces
-(
-    const word& prefix
-) const
-{
-    const pointField& pts = mesh().points();
-
-    const word timeName(Foam::name(mesh().time().value()));
-    const word nameEdges("fluxBased-edges-" + timeName + ".obj");
-    const word nameFaces("fluxBased-faces-" + timeName + ".obj");
-
-
-    // Write OBJ of edge faces to file
-    OBJstream osEdges(mesh().time().path()/nameEdges);
-
-    const auto& edges = mesh().edges();
-    forAll(cornerEdges_, ei)
-    {
-        if (cornerEdges_[ei])
-        {
-            const edge& e = edges[ei];
-            osEdges.write(e, pts);
-        }
-    }
-
-    // Write OBJ of corner faces to file
-    OBJstream osFaces(mesh().time().path()/nameFaces);
-
-    const bitSet& cornerFaces = this->getCornerFaces();
-    const auto& faces = mesh().faces();
-    forAll(cornerFaces, fi)
-    {
-        if (cornerFaces[fi])
-        {
-            const face& f = faces[fi];
-            osFaces.write(f, pts);
-        }
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::cornerDetectionModels::fluxBased::fluxBased
-(
-    const faMesh& mesh,
-    const regionModels::areaSurfaceFilmModels::liquidFilmBase& film,
-    const dictionary& dict
-)
-:
-    cornerDetectionModel(mesh, film, dict),
-    init_(false)
-{
-    read(dict);
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-bool Foam::cornerDetectionModels::fluxBased::detectCorners()
-{
-    if (!init_ || mesh().moving())
-    {
-        // Identify and store corner edges based on face normals
-        cornerEdges_ = identifyCornerEdges();
-        init_ = true;
-    }
-
-
-    // Identify and store corner faces based on edge flux sign
-    this->setCornerFaces(identifyCornerFaces(cornerEdges_));
-
-
-    // Calculate and store corner face angles
-    const bitSet& cornerFaces = this->getCornerFaces();
-    this->setCornerAngles
-    (
-        calcCornerAngles(cornerFaces, cornerEdges_)
-    );
-
-
-    // Write edges and faces as OBJ sets for debug purposes, if need be
-    if (debug && mesh().time().writeTime())
-    {
-        writeEdgesAndFaces();
-    }
-
-    return true;
-}
-
-
-bool Foam::cornerDetectionModels::fluxBased::read(const dictionary& dict)
-{
-    if (!cornerDetectionModel::read(dict))
-    {
-        return false;
-    }
-
-    // Force the re-identification of corner edges/faces
-    init_ = false;
-
-    return true;
-}
-
-
-// ************************************************************************* //

src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/cornerDetectionModels/fluxBased/fluxBased.H

@@ -1,160 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2025 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::cornerDetectionModels::fluxBased
-
-Description
-    Flux-based corner detection model. Marks faces at which the liquid film can
-    separate based on the flux.
-
-    The model identifies sharp edges based on the face normals of the two
-    faces sharing the edge. Then, if the edge is sharp, the flux direction is
-    evaluated to mark the face through which the flux leaves the liquid film.
-    If the edge is sharp and the flux leaves through one of the two faces
-    sharing the edge, the face is marked as a corner face, where the film can
-    separate.
-
-Usage
-    Minimal example in boundary-condition files:
-    \verbatim
-    filmSeparationCoeffs
-    {
-        // Inherited entries
-        ...
-
-        // Optional entries
-        cornerDetectionModel  fluxBased;
-    }
-    \endverbatim
-
-    where the entries mean:
-    \table
-      Property     | Description                     | Type | Reqd | Deflt
-      cornerDetectionModel | Corner detector model   | word | no   | fluxBased
-    \endtable
-
-    The inherited entries are elaborated in:
-      - \link cornerDetectionModel.H \endlink
-
-SourceFiles
-    fluxBased.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef Foam_cornerDetectionModels_fluxBased_H
-#define Foam_cornerDetectionModels_fluxBased_H
-
-#include "cornerDetectionModel.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace cornerDetectionModels
-{
-
-/*---------------------------------------------------------------------------*\
-                          Class fluxBased Declaration
-\*---------------------------------------------------------------------------*/
-
-class fluxBased
-:
-    public cornerDetectionModel
-{
-    // Private Data
-
-        //- Flag to deduce if the object is initialised
-        bool init_;
-
-        //- Identified corner edges
-        bitSet cornerEdges_;
-
-
-    // Private Member Functions
-
-        //- Return Boolean list of identified corner edges
-        bitSet identifyCornerEdges() const;
-
-        //- Return Boolean list of identified corner faces
-        bitSet identifyCornerFaces(const bitSet& cornerEdges) const;
-
-        //- Return the list of corner angles for each edge [rad]
-        scalarList calcCornerAngles
-        (
-            const bitSet& faces,
-            const bitSet& edges
-        ) const;
-
-        // Write edges and faces as OBJ sets for debug purposes
-        void writeEdgesAndFaces(const word& prefix = "geomCorners") const;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("fluxBased");
-
-
-    // Constructors
-
-        //- Construct from components
-        fluxBased
-        (
-            const faMesh& mesh,
-            const regionModels::areaSurfaceFilmModels::liquidFilmBase& film,
-            const dictionary& dict
-        );
-
-
-    // Destructor
-    virtual ~fluxBased() = default;
-
-
-    // Member Functions
-
-    // Evaluation
-
-        //- Detect and store the corner faces
-        virtual bool detectCorners();
-
-
-    // I-O
-
-        //- Read the model dictionary
-        virtual bool read(const dictionary&);
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace cornerDetectionModels
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/cornerDetectionModels/geometryBased/geometryBased.C

@@ -1,586 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2025 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "geometryBased.H"
-#include "processorFaPatch.H"
-#include "unitConversion.H"
-#include "syncTools.H"
-#include "OBJstream.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace cornerDetectionModels
-{
-    defineTypeNameAndDebug(geometryBased, 0);
-    addToRunTimeSelectionTable(cornerDetectionModel, geometryBased, dictionary);
-}
-}
-
-const Foam::Enum
-<
-    Foam::cornerDetectionModels::geometryBased::cornerCurveType
->
-Foam::cornerDetectionModels::geometryBased::cornerCurveTypeNames
-({
-    { cornerCurveType::ANY, "any" },
-    { cornerCurveType::CONCAVE , "concave" },
-    { cornerCurveType::CONVEX , "convex" }
-});
-
-const Foam::Enum
-<
-    Foam::cornerDetectionModels::geometryBased::cornerType
->
-Foam::cornerDetectionModels::geometryBased::cornerTypeNames
-({
-    { cornerType::ALL, "sharpOrRound" },
-    { cornerType::SHARP , "sharp" },
-    { cornerType::ROUND , "round" }
-});
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-Foam::scalar Foam::cornerDetectionModels::geometryBased::curvatureSign
-(
-    const vector& t,
-    const vector& n0,
-    const vector& n1
-) const
-{
-    // t: unit edge tangent
-    // n0: owner face unit normal
-    // n1: neighbour face unit normal
-    scalar curvature = (t & (n0 ^ n1));
-
-    // Orientation: sign of triple product t . (n0 x n1)
-    // Positive => one sense (together with outward normals, treat as "convex");
-    // mapping to convex/concave is finally gated by 'cornerCurveType_'.
-    return sign(curvature);
-}
-
-
-void Foam::cornerDetectionModels::geometryBased::classifyEdges
-(
-    bitSet& sharpEdges,
-    bitSet& roundEdges
-) const
-{
-    // Cache the operand references
-    const areaVectorField& nf = mesh().faceAreaNormals();
-    const edgeList& edges = mesh().edges();
-    const labelUList& own = mesh().edgeOwner();      // own.sz = nEdges
-    const labelUList& nei = mesh().edgeNeighbour();  // nei.sz = nInternalEdges
-    const pointField& pts = mesh().points();
-
-
-    // Convert input-angle parameters from degrees to radians
-    const scalar angSharp = degToRad(angleSharpDeg_);
-    const scalar angRoundMin = degToRad(angleRoundMinDeg_);
-    const scalar angRoundMax = degToRad(angleRoundMaxDeg_);
-
-
-    // Limit to subset of patches if requested
-    bitSet allowedFaces(mesh().nFaces(), true);
-
-
-    // Allocate the resource for the return objects
-    sharpEdges.resize(mesh().nEdges());  // internal + boundary edges
-    sharpEdges.reset();
-
-    roundEdges.resize(mesh().nEdges());
-    roundEdges.reset();
-
-
-    // Internal edges
-    const label nInternalEdges = mesh().nInternalEdges();
-    for (label ei = 0; ei < nInternalEdges; ++ei)
-    {
-        // Do not allow processing of edges shorter than 'minEdgeLength'
-        const edge& e = edges[ei];
-        const scalar le = e.mag(pts);
-        if (le <= max(minEdgeLength_, VSMALL)) continue;
-
-
-        // Do not allow processing of excluded faces
-        const label f0 = own[ei];
-        const label f1 = nei[ei];
-        if (!allowedFaces.test(f0) && !allowedFaces.test(f1)) continue;
-
-
-        // Calculate the dihedral angle and curvature per edge
-        const vector& n0 = nf[f0];
-        const vector& n1 = nf[f1];
-
-        const scalar phi = this->dihedralAngle(n0, n1); // [rad]
-        const scalar kappa = 2.0*Foam::sin(0.5*phi)/max(le, VSMALL); // [1/m]
-        const scalar R = (kappa > VSMALL ? scalar(1)/kappa : GREAT);
-
-        const vector tangent(e.unitVec(pts));
-
-        const scalar sgn = curvatureSign(tangent, n0, n1);
-
-        const bool curvatureType =
-            (cornerCurveType_ == cornerCurveType::ANY)
-         || (cornerCurveType_ == cornerCurveType::CONVEX  && sgn > 0)
-         || (cornerCurveType_ == cornerCurveType::CONCAVE && sgn < 0);
-
-        // Do not allow processing of excluded curvature-type faces
-        if (!curvatureType) continue;
-
-
-        // Sharp: dihedral above threshold
-        if (phi >= angSharp && cornerType_ != cornerType::ROUND)
-        {
-            sharpEdges.set(ei);
-            continue; // do not double-classify as round
-        }
-
-
-        // Round: small-to-moderate angle but small radius (tight fillet)
-        if
-        (
-            phi >= angRoundMin && phi <= angRoundMax
-         && R <= maxRoundRadius_
-         && cornerType_ != cornerType::SHARP
-        )
-        {
-            roundEdges.set(ei);
-        }
-    }
-
-
-    // Optional binary smoothing (edge-neighbour OR)
-
-
-    // Skip the rest of the routine if the simulation is a serial run
-    if (!Pstream::parRun()) return;
-
-    // Boundary edges
-    const faBoundaryMesh& patches = mesh().boundary();
-    for (const faPatch& fap : patches)
-    {
-        const label patchi = fap.index();
-        const label boundaryEdge0 = fap.start();
-
-        const auto& nfp = nf.boundaryField()[patchi];
-
-        if (isA<processorFaPatch>(fap))
-        {
-            forAll(nfp, bEdgei)
-            {
-                const label meshEdgei = boundaryEdge0 + bEdgei;
-
-                // Do not allow processing of edges shorter than 'minEdgeLength'
-                const edge& e = edges[meshEdgei];
-                const scalar le = e.mag(pts);
-                if (le <= max(minEdgeLength_, VSMALL)) continue;
-
-
-                // Do not allow processing of excluded faces
-                const label faceO = own[meshEdgei];
-                if (!allowedFaces.test(faceO)) continue;
-
-
-                // Fetch normal vector of owner and neigh faces
-                const vector& n0 = nf[faceO];
-                const vector& n1 = nfp[bEdgei];
-
-
-                // Calculate the dihedral angle and curvature per edge
-                const scalar phi = this->dihedralAngle(n0, n1); // [rad]
-                const scalar kappa = 2.0*Foam::sin(0.5*phi)/max(le, VSMALL);
-                const scalar R = (kappa > VSMALL ? scalar(1)/kappa : GREAT);
-
-                const vector tangent(e.unitVec(pts));
-
-                const scalar sgn = curvatureSign(tangent, n0, n1);
-
-                const bool curvatureType =
-                    (cornerCurveType_ == cornerCurveType::ANY)
-                 || (cornerCurveType_ == cornerCurveType::CONVEX  && sgn > 0)
-                 || (cornerCurveType_ == cornerCurveType::CONCAVE && sgn < 0);
-
-                // Do not allow processing of excluded curvature-type faces
-                if (!curvatureType) continue;
-
-
-                // Sharp: dihedral above threshold
-                if (phi >= angSharp && cornerType_ != cornerType::ROUND)
-                {
-                    sharpEdges.set(meshEdgei);
-                    continue; // do not double-classify as round
-                }
-
-
-                // Round: small-to-moderate angle but small radius
-                if
-                (
-                    phi >= angRoundMin && phi <= angRoundMax
-                 && R <= maxRoundRadius_
-                 && cornerType_ != cornerType::SHARP
-                )
-                {
-                    roundEdges.set(meshEdgei);
-                }
-            }
-        }
-        else
-        {
-            forAll(nfp, bEdgei)
-            {
-                const label meshEdgei = boundaryEdge0 + bEdgei;
-                const label faceO = own[meshEdgei];
-
-                if (sharpBoundaryEdges_ && allowedFaces.test(faceO))
-                {
-                    sharpEdges.set(meshEdgei);
-                }
-                // Do not allow round edges on physical boundaries
-            }
-        }
-    }
-}
-
-
-void Foam::cornerDetectionModels::geometryBased::edgesToFaces
-(
-    const bitSet& edgeMask,
-    bitSet& faceMask
-) const
-{
-    // Cache the operand references
-    const labelUList& own = mesh().edgeOwner();
-    const labelUList& nei = mesh().edgeNeighbour();
-
-
-    // Allocate the resource for the return objects
-    faceMask.resize(mesh().nFaces());
-    faceMask.reset();
-
-
-    // Internal edges
-    const label nInternalEdges = mesh().nInternalEdges();
-    for (label ei = 0; ei < nInternalEdges; ++ei)
-    {
-        if (edgeMask.test(ei))
-        {
-            // pick the intersecting owner and neighbour faces at the edge
-            faceMask.set(nei[ei]);
-        }
-    }
-
-
-    // Skip the rest of the routine if the simulation is a serial run
-    if (!Pstream::parRun()) return;
-
-
-    // Boundary edges
-    const faBoundaryMesh& patches = mesh().boundary();
-    for (const faPatch& fap : patches)
-    {
-        const label bEdge0 = fap.start();
-        const label nbEdges = fap.size();
-
-        for (label bEdgei = 0; bEdgei < nbEdges; ++bEdgei)
-        {
-            const label meshEdgei = bEdge0 + bEdgei;
-
-            if (edgeMask.test(meshEdgei))
-            {
-                faceMask.set(own[meshEdgei]);
-            }
-        }
-    }
-}
-
-
-Foam::scalarList Foam::cornerDetectionModels::geometryBased::calcCornerAngles
-(
-    const bitSet& faces,
-    const bitSet& edges
-) const
-{
-    // Cache the operand references
-    const areaVectorField& nf = mesh().faceAreaNormals();
-    const labelUList& own = mesh().edgeOwner();
-    const labelUList& nei = mesh().edgeNeighbour();
-
-
-    // Allocate the resource for the return object
-    scalarList cornerFaceAngles(mesh().faces().size(), Zero);
-
-
-    // Internal edges
-    const label nInternalEdges = mesh().nInternalEdges();
-    for (label ei = 0; ei < nInternalEdges; ++ei)
-    {
-        if (!edges[ei]) continue;
-
-        const label faceO = own[ei];
-        const label faceN = nei[ei];
-
-        // If neither adjacent face is flagged as a corner, skip the atan2 work
-        if (!faces[faceO] && !faces[faceN]) continue;
-
-        const scalar ang = this->dihedralAngle(nf[faceO], nf[faceN]);
-
-        if (faces[faceO]) cornerFaceAngles[faceO] = ang;
-        if (faces[faceN]) cornerFaceAngles[faceN] = ang;
-    }
-
-
-    // Skip the rest of the routine if the simulation is a serial run
-    if (!Pstream::parRun()) return cornerFaceAngles;
-
-
-    // Boundary edges
-    const faBoundaryMesh& patches = mesh().boundary();
-    for (const faPatch& fap : patches)
-    {
-        if (!isA<processorFaPatch>(fap)) continue;
-
-        const label patchi = fap.index();
-        const label bEdge0 = fap.start();
-
-        const auto& nfp = nf.boundaryField()[patchi];
-
-        forAll(nfp, bEdgei)
-        {
-            const label meshEdgei = bEdge0 + bEdgei;
-            const label faceO = own[meshEdgei];
-
-            // Only if the mesh edge is a corner and the owner face is a corner
-            if (!edges[meshEdgei] || !faces[faceO]) continue;
-
-            cornerFaceAngles[faceO] =
-                this->dihedralAngle(nf[faceO], nfp[bEdgei]);
-        }
-    }
-
-    return cornerFaceAngles;
-}
-
-
-void Foam::cornerDetectionModels::geometryBased::writeEdgesAndFaces() const
-{
-    // Cache the operand references
-    const auto& edges = mesh().edges();
-    const auto& faces = mesh().faces();
-    const pointField& pts = mesh().points();
-
-    // Generic writer for masked primitives (edge/face)
-    auto writeMasked =
-    [&](
-        const auto& geom,
-        const auto& mask,
-        const char* file
-    )
-    {
-        OBJstream os(mesh().time().path()/file);
-        forAll(mask, i) if (mask[i]) os.write(geom[i], pts);
-    };
-
-    const bool writeSharp =
-        (cornerType_ == cornerType::ALL || cornerType_ == cornerType::SHARP);
-    const bool writeRound =
-        (cornerType_ == cornerType::ALL || cornerType_ == cornerType::ROUND);
-
-    if (writeSharp)
-    {
-        writeMasked(edges, sharpEdges_, "sharp-edges.obj");
-        writeMasked(faces, sharpFaces_, "sharp-faces.obj");
-    }
-
-    if (writeRound)
-    {
-        writeMasked(edges, roundEdges_, "round-edges.obj");
-        writeMasked(faces, roundFaces_, "round-faces.obj");
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::cornerDetectionModels::geometryBased::geometryBased
-(
-    const faMesh& mesh,
-    const regionModels::areaSurfaceFilmModels::liquidFilmBase& film,
-    const dictionary& dict
-)
-:
-    cornerDetectionModel(mesh, film, dict),
-    init_(false)
-{
-    read(dict);
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-bool Foam::cornerDetectionModels::geometryBased::detectCorners()
-{
-    if (!init_ || mesh().moving())
-    {
-        // Identify and store sharp/round edges/faces
-        classifyEdges(sharpEdges_, roundEdges_);
-        edgesToFaces(sharpEdges_, sharpFaces_);
-        edgesToFaces(roundEdges_, roundFaces_);
-
-
-        // Collect all operand edges/faces
-        cornerEdges_ = (sharpEdges_ | roundEdges_);
-        cornerFaces_ = (sharpFaces_ | roundFaces_);
-
-
-        // Pass the operand edges/faces to the film-separation model
-        this->setCornerFaces(cornerFaces_);
-        this->setCornerAngles
-        (
-            calcCornerAngles(cornerFaces_, cornerEdges_)
-        );
-
-
-        init_ = true;
-    }
-
-
-    // Write edges and faces as OBJ sets for debug purposes, if need be
-    if (debug && mesh().time().writeTime())
-    {
-        writeEdgesAndFaces();
-    }
-
-    return true;
-}
-
-
-bool Foam::cornerDetectionModels::geometryBased::read(const dictionary& dict)
-{
-    if (!cornerDetectionModel::read(dict))
-    {
-        return false;
-    }
-
-
-    cornerCurveType_ = cornerCurveTypeNames.getOrDefault
-    (
-        "cornerCurveType",
-        dict,
-        cornerCurveType::ANY
-    );
-
-    cornerType_ = cornerTypeNames.getOrDefault
-    (
-        "cornerType",
-        dict,
-        cornerType::ALL
-    );
-
-    angleSharpDeg_ = dict.getOrDefault<scalar>("angleSharp", 45);
-    angleRoundMinDeg_ = dict.getOrDefault<scalar>("angleRoundMin", 5);
-    angleRoundMaxDeg_ = dict.getOrDefault<scalar>("angleRoundMax", 45);
-    maxRoundRadius_ = dict.getOrDefault<scalar>("maxRoundRadius", 2e-3);
-
-    minEdgeLength_ = dict.getOrDefault<scalar>("minEdgeLength", 0);
-    nSmooth_ = dict.getOrDefault<label>("nSmooth", 0);
-
-    sharpBoundaryEdges_ = dict.getOrDefault<bool>("sharpBoundaryEdges", false);
-
-
-    // Validate the input parameters
-    if (angleSharpDeg_ <= 0 || angleSharpDeg_ >= 180)
-    {
-        FatalIOErrorInFunction(dict)
-            << "angleSharp (" << angleSharpDeg_
-            << " deg) must be in (0, 180)."
-            << exit(FatalIOError);
-    }
-
-    if
-    (
-        angleRoundMinDeg_ < 0 || angleRoundMaxDeg_ > 180
-     || angleRoundMinDeg_ > angleRoundMaxDeg_
-    )
-    {
-        FatalIOErrorInFunction(dict)
-            << "Inconsistent round-angle range: angleRoundMin="
-            << angleRoundMinDeg_ << " deg, angleRoundMax=" << angleRoundMaxDeg_
-            << " deg. Require 0 <= min <= max <= 180."
-            << exit(FatalIOError);
-    }
-
-    if (angleSharpDeg_ <= angleRoundMaxDeg_)
-    {
-        WarningInFunction
-            << "angleSharp (" << angleSharpDeg_
-            << " deg) <= angleRoundMax (" << angleRoundMaxDeg_
-            << " deg): sharp vs round thresholds overlap; "
-            << "classification may be ambiguous."
-            << nl;
-    }
-
-    if (maxRoundRadius_ < 0)
-    {
-        FatalIOErrorInFunction(dict)
-            << "maxRoundRadius must be non-negative."
-            << exit(FatalIOError);
-    }
-
-    if (minEdgeLength_ < 0)
-    {
-        FatalIOErrorInFunction(dict)
-            << "minEdgeLength must be non-negative."
-            << exit(FatalIOError);
-    }
-
-    if (nSmooth_ < 0)
-    {
-        FatalIOErrorInFunction(dict)
-            << "nSmooth must be non-negative."
-            << exit(FatalIOError);
-    }
-
-
-    sharpEdges_.clear();
-    roundEdges_.clear();
-    cornerEdges_.clear();
-    sharpFaces_.clear();
-    roundFaces_.clear();
-    cornerFaces_.clear();
-
-
-    // Force the re-identification of corner edges/faces
-    init_ = false;
-
-    return true;
-}
-
-// ************************************************************************* //

src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/cornerDetectionModels/geometryBased/geometryBased.H

@@ -1,278 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2025 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::cornerDetectionModels::geometryBased
-
-Description
-    Geometry-based corner detection model. Marks faces at which the liquid film
-    can separate based on the geometry.
-
-    The model identifies sharp edges based on the face normals of the two
-    faces sharing the edge. Then, if the edge is sharp, the curvature is
-    evaluated to mark the face through which the flux leaves the liquid film.
-    If the edge is sharp and the curvature is of the specified type, the face
-    is marked as a corner face, where the film can separate.
-
-Usage
-    Minimal example in boundary-condition files:
-    \verbatim
-    filmSeparationCoeffs
-    {
-        // Inherited entries
-        ...
-
-        // Optional entries
-        cornerDetectionModel  geometryBased;
-        cornerCurveType       <word>;
-        cornerType            <word>;
-        angleSharp            <scalar>;   // [deg]
-        angleRoundMin         <scalar>;   // [deg]
-        angleRoundMax         <scalar>;   // [deg]
-        maxRoundRadius        <scalar>;   // [m]
-        minEdgeLength         <scalar>;   // [m]
-        nSmooth               <label>;    // [no. of passes]
-        sharpBoundaryEdges    <bool>;
-    }
-    \endverbatim
-
-    where the entries mean:
-    \table
-      Property     | Description                   | Type | Reqd | Deflt
-      cornerDetectionModel | Corner detector model | word | no   | fluxBased
-      cornerCurveType | Corner-curvature type      | word | no   | any
-      cornerType      | Corner type                | word | no   | sharpOrRound
-      angleSharp      | Sharp-angle limit [deg]    | scalar | no | 45
-      angleRoundMin   | Minimum round-angle limit [deg] | scalar | no | 5
-      angleRoundMax   | Maximum round-angle limit [deg] | scalar | no | 45
-      maxRoundRadius  | Maximum round-radius limit [m] | scalar | no | 2e-3
-      minEdgeLength   | Minimum edge length [m]     | scalar | no | 0
-      nSmooth         | No. of smoothing passes     | label | no | 0
-      sharpBoundaryEdges | Treat boundary edges as sharp | bool | no | false
-    \endtable
-
-    Options for the \c cornerCurve entry:
-    \verbatim
-        any          | Convex or concave corners
-        convex       | Convex corners
-        concave      | Concave corners
-    \endverbatim
-
-    Options for the \c cornerType entry:
-    \verbatim
-        sharp        | Sharp corners
-        round        | Round corners
-        sharpOrRound | Sharp or round corners
-    \endverbatim
-
-    The inherited entries are elaborated in:
-      - \link cornerDetectionModel.H \endlink
-
-SourceFiles
-    geometryBased.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef Foam_cornerDetectionModels_geometryBased_H
-#define Foam_cornerDetectionModels_geometryBased_H
-
-#include "cornerDetectionModel.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace cornerDetectionModels
-{
-
-/*---------------------------------------------------------------------------*\
-                          Class geometryBased Declaration
-\*---------------------------------------------------------------------------*/
-
-class geometryBased
-:
-    public cornerDetectionModel
-{
-    // Private Enumerations
-
-        //- Options for the corner-curvature type
-        enum cornerCurveType : char
-        {
-            ANY = 0,    //!< "Convex or concave corners"
-            CONVEX,     //!< "Convex corners"
-            CONCAVE     //!< "Concave corners"
-        };
-
-        //- Names for cornerCurveType
-        static const Enum<cornerCurveType> cornerCurveTypeNames;
-
-
-        //- Options for the corner type
-        enum cornerType : char
-        {
-            ALL = 0,    //!< "Sharp or round corners"
-            SHARP,      //!< "Sharp corners"
-            ROUND       //!< "Round corners"
-        };
-
-        //- Names for cornerType
-        static const Enum<cornerType> cornerTypeNames;
-
-
-    // Private Data
-
-        //- Corner-curvature type
-        enum cornerCurveType cornerCurveType_;
-
-        //- Corner type
-        enum cornerType cornerType_;
-
-        //- Flag to deduce if the object is initialised
-        bool init_;
-
-        //- Bitset of edges identified as sharp
-        bitSet sharpEdges_;
-
-        //- Bitset of edges identified as round
-        bitSet roundEdges_;
-
-        //- Bitset of edges identified as a combination of sharp and round
-        bitSet cornerEdges_;
-
-        //- Bitset of faces identified as sharp
-        bitSet sharpFaces_;
-
-        //- Bitset of faces identified as round
-        bitSet roundFaces_;
-
-        //- Bitset of faces identified as a combination of sharp and round
-        bitSet cornerFaces_;
-
-        //- Sharp-angle limit
-        scalar angleSharpDeg_;
-
-        //- Minimum round-angle limit
-        scalar angleRoundMinDeg_;
-
-        //- Maximum round-angle limit
-        scalar angleRoundMaxDeg_;
-
-        //- Maximum round-radius limit
-        scalar maxRoundRadius_;
-
-        //- Minimum edge length; ignore edges shorter than this
-        scalar minEdgeLength_;
-
-        //- Number of smoothing passes on the binary edge mask
-        label nSmooth_;
-
-        //- Flag to treat one-face boundary edges as sharp
-        bool sharpBoundaryEdges_;
-
-
-    // Private Member Functions
-
-        //- Return the signed bending sense, sign(+1/-1) of curvature, across
-        //- an edge with respect to the edge tangent
-        scalar curvatureSign
-        (
-            const vector& t,
-            const vector& n0,
-            const vector& n1
-        ) const;
-
-        //- Classify edges into sharp/round according to dihedral angle and
-        //- inferred radius
-        void classifyEdges
-        (
-            bitSet& sharpEdges,
-            bitSet& roundEdges
-        ) const;
-
-        //- Convert an edge mask to a face mask (face is set if any of its
-        //- edges are set)
-        void edgesToFaces
-        (
-            const bitSet& edgeMask,
-            bitSet& faceMask
-        ) const;
-
-        //- Return the list of corner angles [rad] for each edge
-        scalarList calcCornerAngles
-        (
-            const bitSet& faces,
-            const bitSet& edges
-        ) const;
-
-        // Write edges and faces as OBJ sets for debug purposes
-        void writeEdgesAndFaces() const;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("geometryBased");
-
-
-    // Constructors
-
-        //- Construct from components
-        geometryBased
-        (
-            const faMesh& mesh,
-            const regionModels::areaSurfaceFilmModels::liquidFilmBase& film,
-            const dictionary& dict
-        );
-
-
-    // Destructor
-    virtual ~geometryBased() = default;
-
-
-    // Member Functions
-
-    // Evaluation
-
-        //- Detect and store the corner faces
-        virtual bool detectCorners();
-
-
-    // I-O
-
-        //- Read the model dictionary
-        virtual bool read(const dictionary& dict);
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace cornerDetectionModels
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/FriedrichModel/FriedrichModel.C

@@ -26,7 +26,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "FriedrichModel.H"
-#include "cornerDetectionModel.H"
 #include "processorFaPatch.H"
 #include "unitConversion.H"
 #include "addToRunTimeSelectionTable.H"
@@ -54,6 +53,321 @@ FriedrichModel::separationTypeNames
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+bitSet FriedrichModel::calcCornerEdges() const
+{
+    bitSet cornerEdges(mesh().nEdges(), false);
+
+    const areaVectorField& faceCentres = mesh().areaCentres();
+    const areaVectorField& faceNormals = mesh().faceAreaNormals();
+
+    const labelUList& own = mesh().edgeOwner();
+    const labelUList& nbr = mesh().edgeNeighbour();
+
+    // Check if internal face-normal vectors diverge (no separation)
+    // or converge (separation may occur)
+    forAll(nbr, edgei)
+    {
+        const label faceO = own[edgei];
+        const label faceN = nbr[edgei];
+
+        cornerEdges[edgei] = isCornerEdgeSharp
+        (
+            faceCentres[faceO],
+            faceCentres[faceN],
+            faceNormals[faceO],
+            faceNormals[faceN]
+        );
+    }
+
+
+    // Skip the rest of the routine if the simulation is a serial run
+    if (!Pstream::parRun()) return cornerEdges;
+
+    // Check if processor face-normal vectors diverge (no separation)
+    // or converge (separation may occur)
+    const faBoundaryMesh& patches = mesh().boundary();
+
+    for (const faPatch& fap : patches)
+    {
+        if (isA<processorFaPatch>(fap))
+        {
+            const label patchi = fap.index();
+            const auto& edgeFaces = fap.edgeFaces();
+            const label internalEdgei = fap.start();
+
+            const auto& faceCentresp = faceCentres.boundaryField()[patchi];
+            const auto& faceNormalsp = faceNormals.boundaryField()[patchi];
+
+            forAll(faceNormalsp, bndEdgei)
+            {
+                const label faceO = edgeFaces[bndEdgei];
+                const label meshEdgei = internalEdgei + bndEdgei;
+
+                cornerEdges[meshEdgei] = isCornerEdgeSharp
+                (
+                    faceCentres[faceO],
+                    faceCentresp[bndEdgei],
+                    faceNormals[faceO],
+                    faceNormalsp[bndEdgei]
+                );
+            }
+        }
+    }
+
+    return cornerEdges;
+}
+
+
+bool FriedrichModel::isCornerEdgeSharp
+(
+    const vector& faceCentreO,
+    const vector& faceCentreN,
+    const vector& faceNormalO,
+    const vector& faceNormalN
+) const
+{
+    // Calculate the relative position of centres of faces sharing an edge
+    const vector relativePosition(faceCentreN - faceCentreO);
+
+    // Calculate the relative normal of faces sharing an edge
+    const vector relativeNormal(faceNormalN - faceNormalO);
+
+    // Return true if the face normals converge, meaning that the edge is sharp
+    return ((relativeNormal & relativePosition) < -1e-8);
+}
+
+
+scalarList FriedrichModel::calcCornerAngles() const
+{
+    scalarList cornerAngles(mesh().nEdges(), Zero);
+
+    const areaVectorField& faceNormals = mesh().faceAreaNormals();
+
+    const labelUList& own = mesh().edgeOwner();
+    const labelUList& nbr = mesh().edgeNeighbour();
+
+    // Process internal edges
+    forAll(nbr, edgei)
+    {
+        if (!cornerEdges_[edgei]) continue;
+
+        const label faceO = own[edgei];
+        const label faceN = nbr[edgei];
+
+        cornerAngles[edgei] = calcCornerAngle
+        (
+            faceNormals[faceO],
+            faceNormals[faceN]
+        );
+    }
+
+
+    // Skip the rest of the routine if the simulation is a serial run
+    if (!Pstream::parRun()) return cornerAngles;
+
+    // Process processor edges
+    const faBoundaryMesh& patches = mesh().boundary();
+
+    for (const faPatch& fap : patches)
+    {
+        if (isA<processorFaPatch>(fap))
+        {
+            const label patchi = fap.index();
+            const auto& edgeFaces = fap.edgeFaces();
+            const label internalEdgei = fap.start();
+
+            const auto& faceNormalsp = faceNormals.boundaryField()[patchi];
+
+            forAll(faceNormalsp, bndEdgei)
+            {
+                const label faceO = edgeFaces[bndEdgei];
+                const label meshEdgei = internalEdgei + bndEdgei;
+
+                if (!cornerEdges_[meshEdgei]) continue;
+
+                cornerAngles[meshEdgei] = calcCornerAngle
+                (
+                    faceNormals[faceO],
+                    faceNormalsp[bndEdgei]
+                );
+            }
+        }
+    }
+
+    return cornerAngles;
+}
+
+
+scalar FriedrichModel::calcCornerAngle
+(
+    const vector& faceNormalO,
+    const vector& faceNormalN
+) const
+{
+    const scalar magFaceNormal = mag(faceNormalO)*mag(faceNormalN);
+
+    // Avoid any potential exceptions during the cosine calculations
+    if (magFaceNormal < SMALL) return 0;
+
+    scalar cosAngle = (faceNormalO & faceNormalN)/magFaceNormal;
+    cosAngle = clamp(cosAngle, -1, 1);
+
+    return std::acos(cosAngle);
+}
+
+
+bitSet FriedrichModel::calcSeparationFaces() const
+{
+    bitSet separationFaces(mesh().faces().size(), false);
+
+    const edgeScalarField& phis = film().phi2s();
+
+    const labelUList& own = mesh().edgeOwner();
+    const labelUList& nbr = mesh().edgeNeighbour();
+
+    // Process internal faces
+    forAll(nbr, edgei)
+    {
+        if (!cornerEdges_[edgei]) continue;
+
+        const label faceO = own[edgei];
+        const label faceN = nbr[edgei];
+
+        isSeparationFace
+        (
+            separationFaces,
+            phis[edgei],
+            faceO,
+            faceN
+        );
+    }
+
+
+    // Skip the rest of the routine if the simulation is a serial run
+    if (!Pstream::parRun()) return separationFaces;
+
+    // Process processor faces
+    const faBoundaryMesh& patches = mesh().boundary();
+
+    for (const faPatch& fap : patches)
+    {
+        if (isA<processorFaPatch>(fap))
+        {
+            const label patchi = fap.index();
+            const auto& edgeFaces = fap.edgeFaces();
+            const label internalEdgei = fap.start();
+
+            const auto& phisp = phis.boundaryField()[patchi];
+
+            forAll(phisp, bndEdgei)
+            {
+                const label faceO = edgeFaces[bndEdgei];
+                const label meshEdgei(internalEdgei + bndEdgei);
+
+                if (!cornerEdges_[meshEdgei]) continue;
+
+                isSeparationFace
+                (
+                    separationFaces,
+                    phisp[bndEdgei],
+                    faceO
+                );
+            }
+        }
+    }
+
+    return separationFaces;
+}
+
+
+void FriedrichModel::isSeparationFace
+(
+    bitSet& separationFaces,
+    const scalar phiEdge,
+    const label faceO,
+    const label faceN
+) const
+{
+    const scalar tol = 1e-8;
+
+    // Assuming there are no sources/sinks at the edge
+    if (phiEdge > tol)  // From owner to neighbour
+    {
+        separationFaces[faceO] = true;
+    }
+    else if ((phiEdge < -tol) && (faceN != -1))  // From neighbour to owner
+    {
+        separationFaces[faceN] = true;
+    }
+}
+
+
+scalarList FriedrichModel::calcSeparationAngles
+(
+    const bitSet& separationFaces
+) const
+{
+    scalarList separationAngles(mesh().faces().size(), Zero);
+
+    const labelUList& own = mesh().edgeOwner();
+    const labelUList& nbr = mesh().edgeNeighbour();
+
+    // Process internal faces
+    forAll(nbr, edgei)
+    {
+        if (!cornerEdges_[edgei]) continue;
+
+        const label faceO = own[edgei];
+        const label faceN = nbr[edgei];
+
+        if (separationFaces[faceO])
+        {
+            separationAngles[faceO] = cornerAngles_[edgei];
+        }
+
+        if (separationFaces[faceN])
+        {
+            separationAngles[faceN] = cornerAngles_[edgei];
+        }
+    }
+
+
+    // Skip the rest of the routine if the simulation is a serial run
+    if (!Pstream::parRun()) return separationAngles;
+
+    // Process processor faces
+    const edgeScalarField& phis = film().phi2s();
+    const faBoundaryMesh& patches = mesh().boundary();
+
+    for (const faPatch& fap : patches)
+    {
+        if (isA<processorFaPatch>(fap))
+        {
+            const label patchi = fap.index();
+            const auto& edgeFaces = fap.edgeFaces();
+            const label internalEdgei = fap.start();
+
+            const auto& phisp = phis.boundaryField()[patchi];
+
+            forAll(phisp, bndEdgei)
+            {
+                const label faceO = edgeFaces[bndEdgei];
+                const label meshEdgei(internalEdgei + bndEdgei);
+
+                if (!cornerEdges_[meshEdgei]) continue;
+
+                if (separationFaces[faceO])
+                {
+                    separationAngles[faceO] = cornerAngles_[meshEdgei];
+                }
+            }
+        }
+    }
+
+    return separationAngles;
+}
+
+
 tmp<scalarField> FriedrichModel::Fratio() const
 {
     const areaVectorField Up(film().Up());
@@ -64,11 +378,10 @@ tmp<scalarField> FriedrichModel::Fratio() const
     const areaScalarField& sigma = film().sigma();
 
     // Identify the faces where separation may occur
-    const bitSet& separationFaces = cornerDetectorPtr_->getCornerFaces();
+    const bitSet separationFaces(calcSeparationFaces());
 
     // Calculate the corner angles corresponding to the separation faces
-    const scalarList& separationAngles = cornerDetectorPtr_->getCornerAngles();
-
+    const scalarList separationAngles(calcSeparationAngles(separationFaces));
 
     // Initialize the force ratio
     auto tFratio = tmp<scalarField>::New(mesh().faces().size(), Zero);
@@ -118,7 +431,7 @@ tmp<scalarField> FriedrichModel::Fratio() const
         if (isA<processorFaPatch>(fap))
         {
             const label patchi = fap.index();
-            const auto& edgeFaces = fap.edgeFaces();
+            const label internalEdgei = fap.start();
 
             const auto& hp = h.boundaryField()[patchi];
             const auto& Ufp = Uf.boundaryField()[patchi];
@@ -132,18 +445,18 @@ tmp<scalarField> FriedrichModel::Fratio() const
                 // Skip the routine if the face is not a candidate for separation
                 if (!separationFaces[i]) continue;
 
-                const label faceO = edgeFaces[i];
+                const label meshEdgei = internalEdgei + i;
 
                 // Calculate the corner-angle trigonometric values
-                const scalar sinAngle = std::sin(separationAngles[faceO]);
-                const scalar cosAngle = std::cos(separationAngles[faceO]);
+                const scalar sinAngle = std::sin(cornerAngles_[meshEdgei]);
+                const scalar cosAngle = std::cos(cornerAngles_[meshEdgei]);
 
                 // Reynolds number (FLW:Eq. 16)
                 const scalar Re = hp[i]*mag(Ufp[i])*rhop[i]/mup[i];
 
                 // Weber number (FLW:Eq. 17)
-                const vector Urel(Upp[i] - Ufp[i]);
-                const scalar We = hp[i]*rhop_*sqr(mag(Urel))/(2.0*sigmap[i]);
+                const vector Urelp(Upp[i] - Ufp[i]);
+                const scalar We = hp[i]*rhop_*sqr(mag(Urelp))/(2.0*sigmap[i]);
 
                 // Characteristic breakup length (FLW:Eq. 15)
                 const scalar Lb =
@@ -186,12 +499,13 @@ FriedrichModel::FriedrichModel
             separationType::FULL
         )
     ),
-    cornerDetectorPtr_(cornerDetectionModel::New(mesh(), film, dict)),
     rhop_(dict.getScalar("rhop")),
     magG_(mag(film.g().value())),
     C0_(dict.getOrDefault<scalar>("C0", 0.882)),
     C1_(dict.getOrDefault<scalar>("C1", -1.908)),
-    C2_(dict.getOrDefault<scalar>("C2", 1.264))
+    C2_(dict.getOrDefault<scalar>("C2", 1.264)),
+    cornerEdges_(calcCornerEdges()),
+    cornerAngles_(calcCornerAngles())
 {
     if (rhop_ < VSMALL)
     {
@@ -209,18 +523,10 @@ FriedrichModel::FriedrichModel
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-FriedrichModel::~FriedrichModel()
-{}  // cornerDetectionModel was forward declared
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 tmp<scalarField> FriedrichModel::separatedMassRatio() const
 {
-    cornerDetectorPtr_->detectCorners();
-
     tmp<scalarField> tFratio = Fratio();
     const auto& Fratio = tFratio.cref();
 
@@ -270,25 +576,6 @@ tmp<scalarField> FriedrichModel::separatedMassRatio() const
             areaFratio.primitiveFieldRef() = Fratio;
             areaFratio.write();
         }
-
-        {
-            areaScalarField cornerAngles
-            (
-                mesh().newIOobject("cornerAngles"),
-                mesh(),
-                dimensionedScalar(dimless, Zero)
-            );
-
-            const bitSet& cornerFaces = cornerDetectorPtr_->getCornerFaces();
-            const scalarList& angles = cornerDetectorPtr_->getCornerAngles();
-
-            forAll(cornerFaces, i)
-            {
-                if (!cornerFaces[i]) continue;
-                cornerAngles[i] = radToDeg(angles[i]);
-            }
-            cornerAngles.write();
-        }
     }
 
 
@@ -296,16 +583,6 @@ tmp<scalarField> FriedrichModel::separatedMassRatio() const
 }
 
 
-/*
-bool FriedrichModel::read(const dictionary& dict) const
-{
-    // Add the base-class reading later
-    // Read the film separation model dictionary
-
-    return true;
-}
-*/
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace filmSeparationModels

src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/FriedrichModel/FriedrichModel.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2024-2025 OpenCFD Ltd.
+    Copyright (C) 2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -116,14 +116,11 @@ Usage
     filmSeparationCoeffs
     {
         // Mandatory entries
-        model                 Friedrich;
-        rhop                  <scalar>;
+        model               Friedrich;
+        rhop                <scalar>;
 
         // Optional entries
-        separationType        <word>;
-
-        // Inherited entries
-        cornerDetectionModel  <word>;
+        separationType      <word>;
     }
     \endverbatim
 
@@ -133,23 +130,14 @@ Usage
       model        | Model name: Friedrich              | word | yes  | -
       rhop         | Primary-phase density            | scalar | yes  | -
       separationType | Film separation type             | word | no   | full
-      cornerDetectionModel | Corner detector model      | word | no   | flux
     \endtable
 
-    The inherited entries are elaborated in:
-      - \link cornerDetectionModel.H \endlink
-
     Options for the \c separationType entry:
     \verbatim
       full          | Full film separation (Friedrich et al., 2008)
       partial       | Partial film separation (Zhang et al., 2018)
     \endverbatim
 
-    Options for the \c cornerDetectionModel entry:
-    \verbatim
-      flux          | Flux-based corner detection algorithm
-    \endverbatim
-
 SourceFiles
     FriedrichModel.C
 
@@ -164,10 +152,6 @@ SourceFiles
 
 namespace Foam
 {
-
-// Forward Declarations
-class cornerDetectionModel;
-
 namespace filmSeparationModels
 {
 
@@ -197,9 +181,6 @@ class FriedrichModel
         //- Film separation type
         enum separationType separation_;
 
-        //- Corner-detection model
-        mutable autoPtr<cornerDetectionModel> cornerDetectorPtr_;
-
         //- Approximate uniform mass density of primary phase
         scalar rhop_;
 
@@ -215,9 +196,53 @@ class FriedrichModel
         //- Empirical constant for the partial separation model
         scalar C2_;
 
+        //- List of flags identifying sharp-corner edges where separation
+        //- may occur
+        bitSet cornerEdges_;
+
+        //- Corner angles of sharp-corner edges where separation may occur
+        scalarList cornerAngles_;
+
 
     // Private Member Functions
 
+        //- Return Boolean list of identified sharp-corner edges
+        bitSet calcCornerEdges() const;
+
+        //- Return true if the given edge is identified as sharp
+        bool isCornerEdgeSharp
+        (
+            const vector& faceCentreO,
+            const vector& faceCentreN,
+            const vector& faceNormalO,
+            const vector& faceNormalN
+        ) const;
+
+        //- Return the list of sharp-corner angles for each edge
+        scalarList calcCornerAngles() const;
+
+        //- Return the sharp-corner angle for a given edge
+        scalar calcCornerAngle
+        (
+            const vector& faceNormalO,
+            const vector& faceNormalN
+        ) const;
+
+        //- Return Boolean list of identified separation faces
+        bitSet calcSeparationFaces() const;
+
+        //- Return true if the given face is identified as a separation face
+        void isSeparationFace
+        (
+            bitSet& separationFaces,
+            const scalar phiEdge,
+            const label faceO,
+            const label faceN = -1
+        ) const;
+
+        //- Return the list of sharp-corner angles for each face
+        scalarList calcSeparationAngles(const bitSet& separationFaces) const;
+
         //- Return the film-separation force ratio per face
         tmp<scalarField> Fratio() const;
 
@@ -239,7 +264,7 @@ public:
 
 
     // Destructor
-    virtual ~FriedrichModel();
+    virtual ~FriedrichModel() = default;
 
 
     // Member Functions

tutorials/incompressible/simpleFoam/pitzDaily_fused/Allrun

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

tutorials/incompressible/simpleFoam/pitzDaily_fused/system/controlDict

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

tutorials/incompressible/simpleFoam/pitzDaily_fused/system/fvSolution

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

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

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