ENH: Friedrich: refactor model to use cornerDetectionModel framework

- replace inline corner detection logic with cornerDetectionModel integration

src/Pstream/mpi/UPstreamRequest.C

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

src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.C

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

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
-$(kinematic)/injectionModel/filmSeparation/filmSeparationModels/filmSeparationModel/filmSeparationModel.C
-$(kinematic)/injectionModel/filmSeparation/filmSeparationModels/filmSeparationModel/filmSeparationModelNew.C
-$(kinematic)/injectionModel/filmSeparation/filmSeparationModels/OwenRyleyModel/OwenRyleyModel.C
-$(kinematic)/injectionModel/filmSeparation/filmSeparationModels/FriedrichModel/FriedrichModel.C
+/* Film separation models */
+
+filmSeparation=$(kinematic)/injectionModel/filmSeparation
+filmSeparationModels=$(filmSeparation)/filmSeparationModels
+$(filmSeparation)/filmSeparation.C
+$(filmSeparationModels)/filmSeparationModel/filmSeparationModel.C
+$(filmSeparationModels)/filmSeparationModel/filmSeparationModelNew.C
+$(filmSeparationModels)/OwenRyleyModel/OwenRyleyModel.C
+$(filmSeparationModels)/FriedrichModel/FriedrichModel.C
+cornerDetectionModels=$(filmSeparation)/cornerDetectionModels
+$(cornerDetectionModels)/cornerDetectionModel/cornerDetectionModel.C
+$(cornerDetectionModels)/cornerDetectionModel/cornerDetectionModelNew.C
+$(cornerDetectionModels)/fluxBased/fluxBased.C
+$(cornerDetectionModels)/geometryBased/geometryBased.C
 
 $(kinematic)/injectionModel/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/force/contactAngleForces/dynamicContactAngleForce/dynamicContactAngleForce.C

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

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

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

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 cosAngle = std::cos(cornerAngles_[meshEdgei]);
+                const scalar sinAngle = std::sin(separationAngles[faceO]);
+                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 scalar We = hp[i]*rhop_*sqr(mag(Urelp))/(2.0*sigmap[i]);
+                const vector Urel(Upp[i] - Ufp[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)),
-    cornerEdges_(calcCornerEdges()),
-    cornerAngles_(calcCornerAngles())
+    C2_(dict.getOrDefault<scalar>("C2", 1.264))
 {
     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;
-        rhop                <scalar>;
+        model                 Friedrich;
+        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