ENH: Friedrich: refactor model to use cornerDetectionModel framework

- replace inline corner detection logic with cornerDetectionModel integration

applications/test/expressionTemplates/Make/files

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

applications/test/expressionTemplates/Make/options

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

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

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

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

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

src/Allwmake

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

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

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

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

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

src/OpenFOAM/dimensionedTypes/dimensionedType/dimensionedType.H

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

src/OpenFOAM/expressionTemplates/GenericExpression.H

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

src/OpenFOAM/expressionTemplates/boolExpression.H

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

src/OpenFOAM/expressionTemplates/dimensionSetExpression.H

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

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

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

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

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

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

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

src/OpenFOAM/memory/tmp/tmp.H

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

src/OpenFOAM/memory/tmp/tmpI.H

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C

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

src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H

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

src/fused/finiteVolume/fusedGaussGrad.C

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

src/fused/finiteVolume/fusedGaussGrad.H

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

src/fused/finiteVolume/fusedGaussLaplacianScheme.C

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

src/fused/finiteVolume/fusedGaussLaplacianScheme.H

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

src/fused/finiteVolume/fusedGaussLaplacianSchemes.C

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

src/regionFaModels/Make/files

@@ -25,11 +25,20 @@ $(kinematic)/injectionModel/injectionModelList/injectionModelList.C
 $(kinematic)/injectionModel/injectionModel/injectionModel.C
 $(kinematic)/injectionModel/injectionModel/injectionModelNew.C
 
-$(kinematic)/injectionModel/filmSeparation/filmSeparation.C
+/* Film separation models */
-$(kinematic)/injectionModel/filmSeparation/filmSeparationModels/filmSeparationModel/filmSeparationModel.C
+
-$(kinematic)/injectionModel/filmSeparation/filmSeparationModels/filmSeparationModel/filmSeparationModelNew.C
+filmSeparation=$(kinematic)/injectionModel/filmSeparation
-$(kinematic)/injectionModel/filmSeparation/filmSeparationModels/OwenRyleyModel/OwenRyleyModel.C
+filmSeparationModels=$(filmSeparation)/filmSeparationModels
-$(kinematic)/injectionModel/filmSeparation/filmSeparationModels/FriedrichModel/FriedrichModel.C
+$(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/BrunDrippingInjection/BrunDrippingInjection.C
 

src/regionFaModels/Make/options

@@ -11,7 +11,8 @@ 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)/transportModels \
+    -I$(LIB_SRC)/fileFormats/lnInclude
 
 LIB_LIBS = \
     -lfiniteVolume \
@@ -24,4 +25,5 @@ LIB_LIBS = \
     -lthermophysicalProperties \
     -lspecie \
     -lfaOptions \
-    -ldistributionModels
+    -ldistributionModels \
+    -lfileFormats

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

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

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

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

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

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

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

@@ -0,0 +1,278 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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,6 +26,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "FriedrichModel.H"
+#include "cornerDetectionModel.H"
 #include "processorFaPatch.H"
 #include "unitConversion.H"
 #include "addToRunTimeSelectionTable.H"
@@ -53,321 +54,6 @@ 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());
@@ -378,10 +64,11 @@ tmp<scalarField> FriedrichModel::Fratio() const
     const areaScalarField& sigma = film().sigma();
 
     // Identify the faces where separation may occur
-    const bitSet separationFaces(calcSeparationFaces());
+    const bitSet& separationFaces = cornerDetectorPtr_->getCornerFaces();
 
     // Calculate the corner angles corresponding to the separation faces
-    const scalarList separationAngles(calcSeparationAngles(separationFaces));
+    const scalarList& separationAngles = cornerDetectorPtr_->getCornerAngles();
+
 
     // Initialize the force ratio
     auto tFratio = tmp<scalarField>::New(mesh().faces().size(), Zero);
@@ -431,7 +118,7 @@ tmp<scalarField> FriedrichModel::Fratio() const
         if (isA<processorFaPatch>(fap))
         {
             const label patchi = fap.index();
-            const label internalEdgei = fap.start();
+            const auto& edgeFaces = fap.edgeFaces();
 
             const auto& hp = h.boundaryField()[patchi];
             const auto& Ufp = Uf.boundaryField()[patchi];
@@ -445,18 +132,18 @@ tmp<scalarField> FriedrichModel::Fratio() const
                 // Skip the routine if the face is not a candidate for separation
                 if (!separationFaces[i]) continue;
 
-                const label meshEdgei = internalEdgei + i;
+                const label faceO = edgeFaces[i];
 
                 // Calculate the corner-angle trigonometric values
-                const scalar sinAngle = std::sin(cornerAngles_[meshEdgei]);
+                const scalar sinAngle = std::sin(separationAngles[faceO]);
-                const scalar cosAngle = std::cos(cornerAngles_[meshEdgei]);
+                const scalar cosAngle = std::cos(separationAngles[faceO]);
 
                 // Reynolds number (FLW:Eq. 16)
                 const scalar Re = hp[i]*mag(Ufp[i])*rhop[i]/mup[i];
 
                 // Weber number (FLW:Eq. 17)
-                const vector Urelp(Upp[i] - Ufp[i]);
+                const vector Urel(Upp[i] - Ufp[i]);
-                const scalar We = hp[i]*rhop_*sqr(mag(Urelp))/(2.0*sigmap[i]);
+                const scalar We = hp[i]*rhop_*sqr(mag(Urel))/(2.0*sigmap[i]);
 
                 // Characteristic breakup length (FLW:Eq. 15)
                 const scalar Lb =
@@ -499,13 +186,12 @@ 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)
     {
@@ -523,10 +209,18 @@ 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();
 
@@ -576,6 +270,25 @@ 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();
+        }
     }
 
 
@@ -583,6 +296,16 @@ 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 OpenCFD Ltd.
+    Copyright (C) 2024-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -116,11 +116,14 @@ Usage
     filmSeparationCoeffs
     {
         // Mandatory entries
-        model               Friedrich;
+        model                 Friedrich;
-        rhop                <scalar>;
+        rhop                  <scalar>;
 
         // Optional entries
-        separationType      <word>;
+        separationType        <word>;
+
+        // Inherited entries
+        cornerDetectionModel  <word>;
     }
     \endverbatim
 
@@ -130,14 +133,23 @@ 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
 
@@ -152,6 +164,10 @@ SourceFiles
 
 namespace Foam
 {
+
+// Forward Declarations
+class cornerDetectionModel;
+
 namespace filmSeparationModels
 {
 
@@ -181,6 +197,9 @@ class FriedrichModel
         //- Film separation type
         enum separationType separation_;
 
+        //- Corner-detection model
+        mutable autoPtr<cornerDetectionModel> cornerDetectorPtr_;
+
         //- Approximate uniform mass density of primary phase
         scalar rhop_;
 
@@ -196,53 +215,9 @@ 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;
 
@@ -264,7 +239,7 @@ public:
 
 
     // Destructor
-    virtual ~FriedrichModel() = default;
+    virtual ~FriedrichModel();
 
 
     // Member Functions

tutorials/incompressible/simpleFoam/pitzDaily_fused/Allrun

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

tutorials/incompressible/simpleFoam/pitzDaily_fused/system/controlDict

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

tutorials/incompressible/simpleFoam/pitzDaily_fused/system/fvSolution

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

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

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