ENH: radiometerProbes: new function object to monitor radiative heat flux at internal points

The radiometer probes enable measurement of incident radiative heat flux
at arbitrary internal points within the domain, useful for radiation sensor
simulation and heat transfer analysis in participating media applications.

src/Pstream/mpi/UPstreamRequest.C

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

src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.C

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

src/functionObjects/utilities/Make/files

@@ -57,6 +57,7 @@ writeDictionary/writeDictionary.C
 writeObjects/writeObjects.C
 
 thermoCoupleProbes/thermoCoupleProbes.C
+radiometerProbes/radiometerProbes.C
 
 syncObjects/syncObjects.C
 

src/functionObjects/utilities/Make/options

@@ -9,7 +9,9 @@ EXE_INC = \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/ODE/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/transportModels/compressible/lnInclude
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/parallel/distributed/lnInclude
 
 LIB_LIBS = \
     -lfiniteVolume \
@@ -22,4 +24,6 @@ LIB_LIBS = \
     -lsampling \
     -lODE \
     -lfluidThermophysicalModels \
-    -lcompressibleTransportModels
+    -lcompressibleTransportModels \
+    -lradiationModels \
+    -ldistributed

src/functionObjects/utilities/radiometerProbes/radiometerProbes.C

@@ -0,0 +1,260 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "radiometerProbes.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+    defineTypeNameAndDebug(radiometerProbes, 0);
+    addToRunTimeSelectionTable(functionObject, radiometerProbes, dictionary);
+}
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::functionObjects::radiometerProbes::writeFileHeader(Ostream& os)
+{
+    const pointField& locs = probeLocations();
+
+    writeCommented(os, "Probe,Location,Normal");
+
+    os  << nl;
+    for (label i = 0; i < szProbes_; ++i)
+    {
+        const vector& loc = locs[i];
+        const vector& n = n_[i];
+
+        os  << '#' << ' ' << i
+            << ',' << loc.x() << ',' << loc.y() << ',' << loc.z()
+            << ',' << n.x()   << ',' << n.y()   << ',' << n.z()
+            << nl;
+    }
+
+    os  << "# Time";
+    for (int i = 0; i < szProbes_; ++i)
+    {
+        os  << ',' << i;
+    }
+    os  << endl;
+
+    writtenHeader_ = true;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::functionObjects::radiometerProbes::radiometerProbes
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    regionFunctionObject(name, runTime, dict),
+    internalProber(mesh_, dict),
+    writeFile(mesh_, name, typeName, dict),
+    dom_(mesh_.lookupObject<radiation::fvDOM>("radiationProperties")),
+    firstIter_(true)
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * Public Member Functions * * * * * * * * * * * //
+
+bool Foam::functionObjects::radiometerProbes::read(const dictionary& dict)
+{
+    if
+    (
+        !(
+            regionFunctionObject::read(dict)
+         && internalProber::read(dict)
+         && writeFile::read(dict)
+        )
+    )
+    {
+        return false;
+    }
+
+
+    // Skip if the radiation model is inactive
+    if (!dom_.radiation())
+    {
+        WarningInFunction
+            << "The radiation model is inactive."
+            << "Skipping the function object " << type() << ' ' << name()
+            << endl;
+        return false;
+    }
+
+    Log << type() << ':' << name() << ": read" << nl << nl;
+
+
+    // Probe locations are read by 'internalProber'
+    szProbes_ = this->size();
+
+    // If/when fvDOM is updated, the 'read' func is assumed to be executed to
+    // update the fvDOM properties
+    nRay_ = dom_.nRay();
+
+    if (!szProbes_ || !nRay_)
+    {
+        FatalIOErrorInFunction(dict)
+            << "size(probe locations): " << szProbes_ << nl
+            << "size(rays): " << nRay_ << nl
+            << "The input size of probe locations and rays cannot be zero."
+            << exit(FatalIOError);
+    }
+
+    // Read and check size consistency of probe normals with probe locations
+    dict.readEntry("probeNormals", n_);
+    if (n_.size() != szProbes_)
+    {
+        FatalIOErrorInFunction(dict)
+            << "size(probe locations): " << szProbes_ << nl
+            << "size(probe normals): " << n_.size() << nl
+            << "The input size of probe locations and normals must match."
+            << exit(FatalIOError);
+    }
+    n_.normalise();
+
+
+    // Pre-compute and cache inner product of 'n_' and 'dAve_', and 'C_'
+    // This simplification of calculation is valid only if I is non-negative
+    n_dAve_.resize(nRay_);
+    C_.resize(nRay_);
+
+    for (label rayi = 0; rayi < nRay_; ++rayi)
+    {
+        const vector& dAvei = dom_.IRay(rayi).dAve();
+
+        scalarList& n_dAveRay = n_dAve_[rayi];
+        boolList& Cray = C_[rayi];
+
+        n_dAveRay.resize(szProbes_, Zero);
+        Cray.resize(szProbes_, false);
+
+        for (label pi = 0; pi < szProbes_; ++pi)
+        {
+            n_dAveRay[pi] = n_[pi] & dAvei;
+
+            if (n_dAveRay[pi] < 0)  // ray entering the probe
+            {
+                Cray[pi] = true;
+            }
+        }
+    }
+
+    qin_.resize(szProbes_);
+
+
+    if (writeFile::canResetFile())
+    {
+        writeFile::resetFile(typeName);
+    }
+
+    if (writeFile::canWriteHeader())
+    {
+        writeFileHeader(file());
+    }
+
+    return true;
+}
+
+
+bool Foam::functionObjects::radiometerProbes::execute()
+{
+    // Skip if there is no probe to sample, or the radiation model is inactive
+    if (!szProbes_ || !dom_.radiation() || !shouldCalcThisStep())
+    {
+        return false;
+    }
+
+    Log << type() << ' ' << name() << ": execute" << nl << nl;
+
+    qin_ = Zero;  // resized in 'read'
+
+    for (label rayi = 0; rayi < nRay_; ++rayi)
+    {
+        // Radiative intensity field for this ray
+        const volScalarField& I = dom_.IRay(rayi).I();
+
+        // Sample radiative intensity ray at probe locations
+        tmp<scalarField> tIp  = internalProber::sample(I);
+        const scalarField& Ip = tIp.cref();
+
+        const scalarList& n_dAveRay = n_dAve_[rayi];
+        const boolList& Cray = C_[rayi];
+
+        // Add incident radiative heat flux per probe location for each ray
+        for (label pi = 0; pi < szProbes_; ++pi)
+        {
+            if (Cray[pi])
+            {
+                qin_[pi] += Ip[pi]*n_dAveRay[pi];
+            }
+        }
+    }
+
+    return true;
+}
+
+
+bool Foam::functionObjects::radiometerProbes::write()
+{
+    // Skip if there is no probe to sample, or the radiation model is inactive
+    if (!szProbes_ || !dom_.radiation() || !shouldCalcThisStep())
+    {
+        return false;
+    }
+
+    Log << type() << ' ' << name() << ": write" << nl << nl;
+
+    if (UPstream::master())
+    {
+        Ostream& os = file();
+
+        os  << mesh_.time().timeOutputValue();
+        for (label pi = 0; pi < szProbes_; ++pi)
+        {
+            os  << ',' << qin_[pi];
+        }
+        os  << endl;
+    }
+
+    firstIter_ = false;
+
+    return true;
+}
+
+
+// ************************************************************************* //

src/functionObjects/utilities/radiometerProbes/radiometerProbes.H

@@ -0,0 +1,200 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::functionObjects::radiometerProbes
+
+Group
+    grpUtilitiesFunctionObjects
+
+Description
+    Probes the incident radiative heat flux, qin, at arbitrary points within a
+    domain.
+
+Usage
+    Minimal example by using \c system/controlDict.functions:
+    \verbatim
+    radiometer
+    {
+        // Mandatory entries
+        type                radiometerProbes;
+        libs                (utilityFunctionObjects);
+        probeLocations      (<vectorList>);
+        probeNormals        (<vectorList>);
+
+        // Inherited entries
+        ...
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property     | Description                        | Type | Reqd | Deflt
+      type         | Type name: radiometerProbes        | word | yes  | -
+      libs         | Library name: utilityFunctionObjects | word | yes  | -
+      probeLocations | Locations of probes              | vectorList | yes | -
+      probeNormals | Normals of specified probes        | vectorList | yes | -
+    \endtable
+
+    The inherited entries are elaborated in:
+     - \link regionFunctionObject.H \endlink
+     - \link internalProber.H \endlink
+     - \link writeFile.H \endlink
+     - \link fvDOM.H \endlink
+
+Note
+  - The function object can only be used with the \c fvDOM radiation model.
+  - The \c solverFreq input of the \c fvDOM model has superiority over
+    \c executeControl and \c writeControl entries.
+
+SourceFiles
+    radiometerProbes.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_functionObjects_radiometerProbes_H
+#define Foam_functionObjects_radiometerProbes_H
+
+#include "regionFunctionObject.H"
+#include "internalProber.H"
+#include "writeFile.H"
+#include "fvDOM.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class radiometerProbes Declaration
+\*---------------------------------------------------------------------------*/
+
+class radiometerProbes
+:
+    public regionFunctionObject,
+    public internalProber,
+    public writeFile
+{
+    // Private Data
+
+        //- Const reference to the underlying radiation model
+        const radiation::fvDOM& dom_;
+
+        //- Normal vectors of the specified probes
+        vectorField n_;
+
+        //- Pre-computed inner product of probe normals (n_) and average
+        //- solid-angle direction (dAve) per radiative intensity ray
+        List<scalarList> n_dAve_;
+
+        //- Directional selection coefficient for radiative intensity rays
+        //  false: ray entering the probe
+        //  true: ray leaving the probe
+        List<boolList> C_;
+
+        //- Incident radiative heat flux per probe location
+        scalarField qin_;
+
+        //- Number of radiative intensity rays
+        label nRay_;
+
+        //- Number of probe locations/normals
+        label szProbes_;
+
+        //- Flag to identify whether the iteration is the first iteration
+        //  Resets with a restarted simulation
+        bool firstIter_;
+
+
+    // Private Member Functions
+
+        //- Write file-header information into the output file
+        virtual void writeFileHeader(Ostream& os);
+
+        //- Return the flag to decide if radiation-model calculations are
+        //- performed, so that function object calculations can proceed
+        bool shouldCalcThisStep() const
+        {
+            return
+                firstIter_
+             || (mesh_.time().timeIndex() % dom_.solverFreq() == 0);
+        }
+
+
+public:
+
+    //- Runtime type information
+    TypeName("radiometerProbes");
+
+
+    // Generated Methods
+
+        //- No copy construct
+        radiometerProbes(const radiometerProbes&) = delete;
+
+        //- No copy assignment
+        void operator=(const radiometerProbes&) = delete;
+
+
+    // Constructors
+
+        //- Construct from name, Time and dictionary
+        radiometerProbes
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~radiometerProbes() = default;
+
+
+    // Public Member Functions
+
+        //- Read the function object settings
+        virtual bool read(const dictionary&);
+
+        //- Execute the function object
+        virtual bool execute();
+
+        //- Write to data files/fields and to streams
+        virtual bool write();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/functionObjects/utilities/thermoCoupleProbes/thermoCoupleProbes.C

@@ -76,7 +76,7 @@ Foam::functionObjects::thermoCoupleProbes::thermoCoupleProbes
     }
     else
     {
-        Ttc_ = probes::sample(thermo_.T());
+        Ttc_ = prober().sample(thermo_.T());
     }
 
     // Note: can only create the solver once all samples have been found
@@ -89,7 +89,7 @@ Foam::functionObjects::thermoCoupleProbes::thermoCoupleProbes
 
 Foam::label Foam::functionObjects::thermoCoupleProbes::nEqns() const
 {
-    return this->size();
+    return prober().size();
 }
 
 
@@ -108,19 +108,21 @@ void Foam::functionObjects::thermoCoupleProbes::derivatives
     scalarField Cpc(y.size(), Zero);
     scalarField kappac(y.size(), Zero);
 
+    const auto& p = prober();
+
     if (radiationFieldName_ != "none")
     {
-        G = sample(mesh_.lookupObject<volScalarField>(radiationFieldName_));
+        G = p.sample(mesh_.lookupObject<volScalarField>(radiationFieldName_));
     }
 
-    Tc = probes::sample(thermo_.T());
+    Tc = p.sample(thermo_.T());
 
-    Uc = mag(this->sample(mesh_.lookupObject<volVectorField>(UName_)));
+    Uc = mag(p.sample(mesh_.lookupObject<volVectorField>(UName_)));
 
-    rhoc = this->sample(thermo_.rho()());
-    kappac = this->sample(thermo_.kappa()());
-    muc = this->sample(thermo_.mu()());
-    Cpc = this->sample(thermo_.Cp()());
+    rhoc = p.sample(thermo_.rho()());
+    kappac = p.sample(thermo_.kappa()());
+    muc = p.sample(thermo_.mu()());
+    Cpc = p.sample(thermo_.Cp()());
 
     scalarField Re(rhoc*Uc*d_/muc);
     scalarField Pr(Cpc*muc/kappac);
@@ -163,7 +165,7 @@ void Foam::functionObjects::thermoCoupleProbes::jacobian
 
 bool Foam::functionObjects::thermoCoupleProbes::write()
 {
-    if (!pointField::empty())
+    if (!prober().empty())
     {
         (void) prepare(ACTION_WRITE);
 
@@ -182,7 +184,7 @@ bool Foam::functionObjects::thermoCoupleProbes::write()
 
 bool Foam::functionObjects::thermoCoupleProbes::execute()
 {
-    if (!pointField::empty())
+    if (!prober().empty())
     {
         scalar dt = mesh_.time().deltaTValue();
         scalar t = mesh_.time().value();

src/functionObjects/utilities/thermoCoupleProbes/thermoCoupleProbesTemplates.C

@@ -42,7 +42,8 @@ void Foam::functionObjects::thermoCoupleProbes::writeValues
 
         os  << setw(w) << timeValue;
 
-        forAll(*this, probei)
+        const pointField& probes = prober().probeLocations();
+        forAll(probes, probei)
         {
             // if (includeOutOfBounds_ || processor_[probei] != -1)
             {

src/regionFaModels/Make/files

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

src/regionFaModels/Make/options

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

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

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020-2025 OpenCFD Ltd.
+    Copyright (C) 2020-2023 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,7 +59,6 @@ dynamicContactAngleForce::dynamicContactAngleForce
 )
 :
     contactAngleForce(typeName, film, dict),
-    thetaPtr_(nullptr),
     U_vs_thetaPtr_
     (
         Function1<scalar>::NewIfPresent
@@ -81,56 +80,46 @@ dynamicContactAngleForce::dynamicContactAngleForce
         )
     ),
     rndGen_(label(0)),
-    distributionPtr_(nullptr)
+    distribution_
+    (
+        distributionModel::New
+        (
+            coeffDict_.subDict("distribution"),
+            rndGen_
+        )
+    )
 {
     if (U_vs_thetaPtr_ && T_vs_thetaPtr_)
     {
         FatalIOErrorInFunction(dict)
-            << "Only one of Utheta or Ttheta should be provided; "
-            << "both inputs cannot be used together."
+            << "Entries Utheta and Ttheta could not 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
 {
-    areaScalarField& theta = thetaPtr_.ref();
+    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();
     scalarField& thetai = theta.ref();
 
+
     if (U_vs_thetaPtr_)
     {
         // Initialize with the function of film speed
@@ -147,16 +136,13 @@ tmp<areaScalarField> dynamicContactAngleForce::theta() const
         thetai = T_vs_thetaPtr_->value(T());
     }
 
-    if (distributionPtr_)
+    // Add the stochastic perturbation
+    forAll(thetai, facei)
     {
-        // Add the stochastic perturbation
-        forAll(thetai, facei)
-        {
-            thetai[facei] += distributionPtr_->sample();
-        }
+        thetai[facei] += distribution_->sample();
     }
 
-    return thetaPtr_();
+    return ttheta;
 }
 
 

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

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020-2025 OpenCFD Ltd.
+    Copyright (C) 2020-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -112,9 +112,6 @@ 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_;
 
@@ -124,8 +121,17 @@ class dynamicContactAngleForce
         //- Random number generator
         Random rndGen_;
 
-        //- Stochastic perturbation model for contact angle
-        autoPtr<distributionModel> distributionPtr_;
+        //- 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;
 
 
 protected:
@@ -142,7 +148,7 @@ public:
 
     // Constructors
 
-        //- Construct from surface film model and dictionary
+        //- Construct from surface film model
         dynamicContactAngleForce
         (
             liquidFilmBase& film,
@@ -151,7 +157,7 @@ public:
 
 
     //- Destructor
-    virtual ~dynamicContactAngleForce();
+    virtual ~dynamicContactAngleForce() = default;
 };
 
 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

src/sampling/Make/files

@@ -1,3 +1,6 @@
+probes/probers/prober.C
+probes/probers/internalProber/internalProber.C
+probes/probers/patchProber/patchProber.C
 probes/probes.C
 probes/patchProbes.C
 

src/sampling/probes/Probes/Probes.C

@@ -0,0 +1,453 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2015-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 "Probes.H"
+#include "IOmanip.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "SpanStream.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Prober>
+void Foam::Probes<Prober>::createProbeFiles(const wordList& fieldNames)
+{
+    // Open new output streams
+
+    bool needsNewFiles = false;
+    for (const word& fieldName : fieldNames)
+    {
+        if (!probeFilePtrs_.found(fieldName))
+        {
+            needsNewFiles = true;
+            break;
+        }
+    }
+
+    if (needsNewFiles && Pstream::master())
+    {
+        DebugInfo
+            << "Probing fields: " << fieldNames << nl
+            << "Probing locations: " << prober_.probeLocations() << nl
+            << endl;
+
+        // Put in undecomposed case
+        // (Note: gives problems for distributed data running)
+
+        fileName probeDir
+        (
+            mesh_.time().globalPath()
+          / functionObject::outputPrefix
+          / name()/mesh_.regionName()
+            // Use startTime as the instance for output files
+          / mesh_.time().timeName(mesh_.time().startTime().value())
+        );
+        probeDir.clean();  // Remove unneeded ".."
+
+        // Create directory if needed
+        Foam::mkDir(probeDir);
+
+        for (const word& fieldName : fieldNames)
+        {
+            if (probeFilePtrs_.found(fieldName))
+            {
+                // Safety
+                continue;
+            }
+
+            auto osPtr = autoPtr<OFstream>::New(probeDir/fieldName);
+            auto& os = *osPtr;
+
+            if (!os.good())
+            {
+                FatalErrorInFunction
+                    << "Cannot open probe output file: " << os.name() << nl
+                    << exit(FatalError);
+            }
+
+            probeFilePtrs_.insert(fieldName, osPtr);
+
+            DebugInfo<< "open probe stream: " << os.name() << endl;
+
+            const unsigned int width(IOstream::defaultPrecision() + 7);
+            os.setf(std::ios_base::left);
+
+            const pointField& probeLocs = prober_.probeLocations();
+            const labelList& processors = prober_.processors();
+            const labelList& patchIDList = prober_.patchIDList();
+            const pointField& oldPoints = prober_.oldPoints();
+
+            forAll(probeLocs, probei)
+            {
+                os  << "# Probe " << probei << ' ' << probeLocs[probei];
+
+                if (processors[probei] == -1)
+                {
+                    os  << "  # Not Found";
+                }
+                else if (probei < patchIDList.size())
+                {
+                    const label patchi = patchIDList[probei];
+                    if (patchi != -1)
+                    {
+                        const polyBoundaryMesh& bm = mesh_.boundaryMesh();
+                        if
+                        (
+                            patchi < bm.nNonProcessor()
+                         || processors[probei] == Pstream::myProcNo()
+                        )
+                        {
+                            os  << " at patch " << bm[patchi].name();
+                        }
+                        os  << " with a distance of "
+                            << mag(probeLocs[probei]-oldPoints[probei])
+                            << " m to the original point "
+                            << oldPoints[probei];
+                    }
+                }
+
+                os  << nl;
+            }
+
+            os  << setw(width) << "# Time";
+
+            forAll(probeLocs, probei)
+            {
+                if (prober_.includeOutOfBounds() || processors[probei] != -1)
+                {
+                    os  << ' ' << setw(width) << probei;
+                }
+            }
+            os  << endl;
+        }
+    }
+}
+
+
+template<class Prober>
+template<class Type>
+void Foam::Probes<Prober>::writeValues
+(
+    const word& fieldName,
+    const Field<Type>& values,
+    const scalar timeValue
+)
+{
+    if (Pstream::master())
+    {
+        const unsigned int width(IOstream::defaultPrecision() + 7);
+        OFstream& os = *probeFilePtrs_[fieldName];
+
+        os  << setw(width) << timeValue;
+
+        OCharStream buf;
+
+        const bool includeOutOfBounds = prober_.includeOutOfBounds();
+        const labelList& procs = prober_.processors();
+        forAll(values, probei)
+        {
+            if (includeOutOfBounds || procs[probei] != -1)
+            {
+                buf.rewind();
+                buf << values[probei];
+                os  << ' ' << setw(width) << buf.str().data();
+            }
+        }
+        os  << endl;
+    }
+}
+
+
+template<class Prober>
+template<class GeoField>
+void Foam::Probes<Prober>::performAction
+(
+    const fieldGroup<GeoField>& fieldNames,
+    unsigned request
+)
+{
+    for (const word& fieldName : fieldNames)
+    {
+        tmp<GeoField> tfield = getOrLoadField<GeoField>(fieldName);
+
+        if (tfield)
+        {
+            const auto& field = tfield();
+            const scalar timeValue = field.time().timeOutputValue();
+
+            Field<typename GeoField::value_type> values(prober_.sample(field));
+
+            this->storeResults(fieldName, values);
+            if (request & ACTION_WRITE)
+            {
+                this->writeValues(fieldName, values, timeValue);
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<class Prober>
+Foam::label Foam::Probes<Prober>::prepare(unsigned request)
+{
+    // Prefilter on selection
+    HashTable<wordHashSet> selected =
+    (
+        loadFromFiles_
+      ? IOobjectList(mesh_, mesh_.time().timeName()).classes(fieldSelection_)
+      : mesh_.classes(fieldSelection_)
+    );
+
+    // Classify and count fields
+    label nFields = 0;
+    do
+    {
+        #undef  doLocalCode
+        #define doLocalCode(InputType, Target)                                \
+        {                                                                     \
+            Target.clear();  /* Remove old values */                          \
+            const auto iter = selected.cfind(InputType::typeName);            \
+            if (iter.good())                                                  \
+            {                                                                 \
+                /* Add new (current) values */                                \
+                Target.append(iter.val().sortedToc());                        \
+                nFields += Target.size();                                     \
+            }                                                                 \
+        }
+
+        doLocalCode(volScalarField, scalarFields_);
+        doLocalCode(volVectorField, vectorFields_);
+        doLocalCode(volSphericalTensorField, sphericalTensorFields_);
+        doLocalCode(volSymmTensorField, symmTensorFields_);
+        doLocalCode(volTensorField, tensorFields_);
+
+        doLocalCode(surfaceScalarField, surfaceScalarFields_);
+        doLocalCode(surfaceVectorField, surfaceVectorFields_);
+        doLocalCode(surfaceSphericalTensorField, surfaceSphericalTensorFields_);
+        doLocalCode(surfaceSymmTensorField, surfaceSymmTensorFields_);
+        doLocalCode(surfaceTensorField, surfaceTensorFields_);
+        #undef doLocalCode
+    }
+    while (false);
+
+
+    // Adjust file streams
+    if (Pstream::master())
+    {
+        wordHashSet currentFields(2*nFields);
+        currentFields.insert(scalarFields_);
+        currentFields.insert(vectorFields_);
+        currentFields.insert(sphericalTensorFields_);
+        currentFields.insert(symmTensorFields_);
+        currentFields.insert(tensorFields_);
+
+        currentFields.insert(surfaceScalarFields_);
+        currentFields.insert(surfaceVectorFields_);
+        currentFields.insert(surfaceSphericalTensorFields_);
+        currentFields.insert(surfaceSymmTensorFields_);
+        currentFields.insert(surfaceTensorFields_);
+
+        DebugInfo
+            << "Probing fields: " << currentFields << nl
+            << "Probing locations: " << prober_.probeLocations() << nl
+            << endl;
+
+        // Close streams for fields that no longer exist
+        forAllIters(probeFilePtrs_, iter)
+        {
+            if (!currentFields.erase(iter.key()))
+            {
+                DebugInfo<< "close probe stream: " << iter()->name() << endl;
+
+                probeFilePtrs_.remove(iter);
+            }
+        }
+
+        if ((request & ACTION_WRITE) && !currentFields.empty())
+        {
+            createProbeFiles(currentFields.sortedToc());
+        }
+    }
+
+    return nFields;
+}
+
+
+template<class Prober>
+template<class GeoField>
+Foam::tmp<GeoField>
+Foam::Probes<Prober>::getOrLoadField(const word& fieldName) const
+{
+    tmp<GeoField> tfield;
+
+    if (loadFromFiles_)
+    {
+        tfield.emplace
+        (
+            IOobject
+            (
+                fieldName,
+                mesh_.time().timeName(),
+                mesh_.thisDb(),
+                IOobjectOption::MUST_READ,
+                IOobjectOption::NO_WRITE,
+                IOobjectOption::NO_REGISTER
+            ),
+            mesh_
+        );
+    }
+    else
+    {
+        tfield.cref(mesh_.cfindObject<GeoField>(fieldName));
+    }
+
+    return tfield;
+}
+
+
+template<class Prober>
+template<class Type>
+void Foam::Probes<Prober>::storeResults
+(
+    const word& fieldName,
+    const Field<Type>& values
+)
+{
+    const MinMax<Type> limits(values);
+    const Type avgVal = average(values);
+
+    this->setResult("average(" + fieldName + ")", avgVal);
+    this->setResult("min(" + fieldName + ")", limits.min());
+    this->setResult("max(" + fieldName + ")", limits.max());
+    this->setResult("size(" + fieldName + ")", values.size());
+
+    if (verbose_)
+    {
+        Info<< name() << " : " << fieldName << nl
+            << "    avg: " << avgVal << nl
+            << "    min: " << limits.min() << nl
+            << "    max: " << limits.max() << nl << nl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Prober>
+Foam::Probes<Prober>::Probes
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict,
+    const bool loadFromFiles,
+    const bool readFields
+)
+:
+    functionObjects::fvMeshFunctionObject(name, runTime, dict),
+    prober_(mesh_, dict),
+    loadFromFiles_(loadFromFiles),
+    onExecute_(false),
+    fieldSelection_()
+{
+    if (readFields)
+    {
+        read(dict);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Prober>
+bool Foam::Probes<Prober>::verbose(const bool on) noexcept
+{
+    bool old(verbose_);
+    verbose_ = on;
+    return old;
+}
+
+
+template<class Prober>
+bool Foam::Probes<Prober>::read(const dictionary& dict)
+{
+    dict.readEntry("fields", fieldSelection_);
+
+    verbose_ = dict.getOrDefault("verbose", false);
+    onExecute_ = dict.getOrDefault("sampleOnExecute", false);
+
+    // Close old (unused) streams
+    prepare(ACTION_NONE);
+
+    return true;
+}
+
+
+template<class Prober>
+bool Foam::Probes<Prober>::performAction(unsigned request)
+{
+    if (!prober_.empty() && request && prepare(request))
+    {
+        performAction(scalarFields_, request);
+        performAction(vectorFields_, request);
+        performAction(sphericalTensorFields_, request);
+        performAction(symmTensorFields_, request);
+        performAction(tensorFields_, request);
+
+        performAction(surfaceScalarFields_, request);
+        performAction(surfaceVectorFields_, request);
+        performAction(surfaceSphericalTensorFields_, request);
+        performAction(surfaceSymmTensorFields_, request);
+        performAction(surfaceTensorFields_, request);
+    }
+    return true;
+}
+
+
+template<class Prober>
+bool Foam::Probes<Prober>::execute()
+{
+    if (onExecute_)
+    {
+        return performAction(ACTION_ALL & ~ACTION_WRITE);
+    }
+
+    return true;
+}
+
+
+template<class Prober>
+bool Foam::Probes<Prober>::write()
+{
+    return performAction(ACTION_ALL);
+}
+
+
+// ************************************************************************* //

src/sampling/probes/Probes/Probes.H

@@ -0,0 +1,238 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2016-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::Probes
+
+Description
+    Base class for sampling fields at specified locations and writing to file.
+    The locations are specified and determined in the derived class. The
+    sampling is done using the specified point prober class.
+
+SourceFiles
+    Probes.C
+    ProbesTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_Probes_H
+#define Foam_Probes_H
+
+#include "fvMeshFunctionObject.H"
+#include "polyMesh.H"
+#include "HashPtrTable.H"
+#include "OFstream.H"
+#include "volFieldsFwd.H"
+#include "surfaceFieldsFwd.H"
+#include "prober.H"
+#include "IOobjectList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class Probes Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Prober>
+class Probes
+:
+    public functionObjects::fvMeshFunctionObject
+{
+protected:
+
+    // Protected Data
+
+        //- The specified point prober
+        Prober prober_;
+
+
+    // Protected Classes
+
+        //- Grouping of field names by GeometricField type
+        template<class GeoField>
+        struct fieldGroup : public DynamicList<word> {};
+
+
+    // Data Types
+
+        //- Local control for sampling actions
+        enum sampleActionType : unsigned
+        {
+            ACTION_NONE  = 0,
+            ACTION_WRITE = 0x1,
+            ACTION_STORE = 0x2,
+            ACTION_ALL = 0xF
+        };
+
+
+    // Protected Data
+
+        //- Load fields from files (not from objectRegistry)
+        bool loadFromFiles_;
+
+        //- Output verbosity
+        bool verbose_;
+
+        //- Perform sample actions on execute as well
+        bool onExecute_;
+
+        //- Requested names of fields to probe
+        wordRes fieldSelection_;
+
+
+    // Calculated
+
+        //- Current list of field names selected for sampling
+        DynamicList<word> selectedFieldNames_;
+
+        //- Categorized scalar/vector/tensor volume fields
+        fieldGroup<volScalarField> scalarFields_;
+        fieldGroup<volVectorField> vectorFields_;
+        fieldGroup<volSphericalTensorField> sphericalTensorFields_;
+        fieldGroup<volSymmTensorField> symmTensorFields_;
+        fieldGroup<volTensorField> tensorFields_;
+
+        //- Categorized scalar/vector/tensor surface fields
+        fieldGroup<surfaceScalarField> surfaceScalarFields_;
+        fieldGroup<surfaceVectorField> surfaceVectorFields_;
+        fieldGroup<surfaceSphericalTensorField> surfaceSphericalTensorFields_;
+        fieldGroup<surfaceSymmTensorField> surfaceSymmTensorFields_;
+        fieldGroup<surfaceTensorField> surfaceTensorFields_;
+
+        //- Current open files (non-empty on master only)
+        HashPtrTable<OFstream> probeFilePtrs_;
+
+
+    // Protected Member Functions
+
+        //- Classify field types, close/open file streams
+        //  \return number of fields to sample
+        label prepare(unsigned request);
+
+        //- Get from registry or load from disk
+        template<class GeoField>
+        tmp<GeoField> getOrLoadField(const word& fieldName) const;
+
+        //- Store results: min/max/average/size
+        template<class Type>
+        void storeResults(const word& fieldName, const Field<Type>& values);
+
+
+private:
+
+    // Private Member Functions
+
+        //- Create new streams as required
+        void createProbeFiles(const wordList& fieldNames);
+
+        //- Write field values
+        template<class Type>
+        void writeValues
+        (
+            const word& fieldName,
+            const Field<Type>& values,
+            const scalar timeValue
+        );
+
+        //- Sample and store/write all applicable sampled fields
+        template<class GeoField>
+        void performAction
+        (
+            const fieldGroup<GeoField>& fieldNames,  /* must be sorted */
+            unsigned request
+        );
+
+        //- Perform sampling action with store/write
+        bool performAction(unsigned request);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from Time and dictionary
+        Probes
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict,
+            const bool loadFromFiles = false,
+            const bool readFields = true
+        );
+
+
+    //- Destructor
+    virtual ~Probes() = default;
+
+
+    // Member Functions
+
+        //- Enable/disable verbose output
+        //  \return old value
+        bool verbose(const bool on) noexcept;
+
+        //- Return names of fields to probe
+        virtual const wordRes& fieldNames() const noexcept
+        {
+            return fieldSelection_;
+        }
+
+        //- Return const reference to the point prober
+        const Prober& prober() const noexcept { return prober_; }
+
+        //- Read the settings from the dictionary
+        virtual bool read(const dictionary&);
+
+        //- Sample and store result if the sampleOnExecute is enabled.
+        virtual bool execute();
+
+        //- Sample and write
+        virtual bool write();
+
+        //- Update for changes of mesh due to readUpdate
+        virtual void readUpdate(const polyMesh::readUpdateState state)
+        {}
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "Probes.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/sampling/probes/patchProbes.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2022 OpenCFD Ltd.
+    Copyright (C) 2016-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,11 +27,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "patchProbes.H"
-#include "volFields.H"
-#include "IOmanip.H"
-#include "mappedPatchBase.H"
-#include "treeBoundBox.H"
-#include "treeDataFace.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -39,7 +34,6 @@ License
 namespace Foam
 {
     defineTypeNameAndDebug(patchProbes, 0);
-
     addToRunTimeSelectionTable
     (
         functionObject,
@@ -48,179 +42,6 @@ namespace Foam
     );
 }
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-void Foam::patchProbes::findElements(const fvMesh& mesh)
-{
-    (void)mesh.tetBasePtIs();
-
-    const polyBoundaryMesh& bm = mesh.boundaryMesh();
-
-    // All the info for nearest. Construct to miss
-    List<mappedPatchBase::nearInfo> nearest(this->size());
-
-    const labelList patchIDs(bm.patchSet(patchNames_).sortedToc());
-
-    label nFaces = 0;
-    forAll(patchIDs, i)
-    {
-        nFaces += bm[patchIDs[i]].size();
-    }
-
-    if (nFaces > 0)
-    {
-        // Collect mesh faces and bounding box
-        labelList bndFaces(nFaces);
-        treeBoundBox overallBb;
-
-        nFaces = 0;
-        forAll(patchIDs, i)
-        {
-            const polyPatch& pp = bm[patchIDs[i]];
-            forAll(pp, i)
-            {
-                bndFaces[nFaces++] = pp.start()+i;
-                const face& f = pp[i];
-
-                // Without reduction.
-                overallBb.add(pp.points(), f);
-            }
-        }
-
-        Random rndGen(123456);
-        overallBb.inflate(rndGen, 1e-4, ROOTVSMALL);
-
-
-        const indexedOctree<treeDataFace> boundaryTree
-        (
-            treeDataFace(mesh, bndFaces),  // patch faces only
-
-            overallBb,      // overall search domain
-            8,              // maxLevel
-            10,             // leafsize
-            3.0             // duplicity
-        );
-
-        forAll(probeLocations(), probei)
-        {
-            const auto& treeData = boundaryTree.shapes();
-            const point sample = probeLocations()[probei];
-
-            pointIndexHit info = boundaryTree.findNearest
-            (
-                sample,
-                Foam::sqr(boundaryTree.bb().mag())
-            );
-
-            if (!info.hit())
-            {
-                info = boundaryTree.findNearest(sample, Foam::sqr(GREAT));
-            }
-
-            const label facei = treeData.objectIndex(info.index());
-
-            const label patchi = bm.whichPatch(facei);
-
-            if (isA<emptyPolyPatch>(bm[patchi]))
-            {
-                WarningInFunction
-                    << " The sample point: " << sample
-                    << " belongs to " << patchi
-                    << " which is an empty patch. This is not permitted. "
-                    << " This sample will not be included "
-                    << endl;
-            }
-            else if (info.hit())
-            {
-                // Note: do we store the face centre or the actual nearest?
-                // We interpolate using the faceI only though (no
-                // interpolation) so it does not actually matter much, just for
-                // the location written to the header.
-
-                //const point& facePt = mesh.faceCentres()[faceI];
-                const point& facePt = info.point();
-
-                mappedPatchBase::nearInfo sampleInfo;
-
-                sampleInfo.first() = pointIndexHit(true, facePt, facei);
-
-                sampleInfo.second().first() = facePt.distSqr(sample);
-                sampleInfo.second().second() = Pstream::myProcNo();
-
-                nearest[probei] = sampleInfo;
-            }
-        }
-    }
-
-
-    // Find nearest - globally consistent
-    Pstream::listCombineReduce(nearest, mappedPatchBase::nearestEqOp());
-
-    oldPoints_.resize(this->size());
-
-    // Update actual probe locations and store old ones
-    forAll(nearest, samplei)
-    {
-        oldPoints_[samplei] = operator[](samplei);
-        operator[](samplei) = nearest[samplei].first().point();
-    }
-
-    if (debug)
-    {
-        InfoInFunction << nl;
-        forAll(nearest, samplei)
-        {
-            label proci = nearest[samplei].second().second();
-            label locali = nearest[samplei].first().index();
-
-            Info<< "    " << samplei << " coord:"<< operator[](samplei)
-                << " found on processor:" << proci
-                << " in local face:" << locali
-                << " with location:" << nearest[samplei].first().point()
-                << endl;
-        }
-    }
-
-    // Extract any local faces to sample:
-    // - operator[] : actual point to sample (=nearest point on patch)
-    // - oldPoints_ : original provided point (might be anywhere in the mesh)
-    // - elementList_   : cells, not used
-    // - faceList_      : faces (now patch faces)
-    // - patchIDList_   : patch corresponding to faceList
-    // - processor_     : processor
-    elementList_.resize_nocopy(nearest.size());
-    elementList_ = -1;
-
-    faceList_.resize_nocopy(nearest.size());
-    faceList_ = -1;
-
-    processor_.resize_nocopy(nearest.size());
-    processor_ = -1;
-
-    patchIDList_.resize_nocopy(nearest.size());
-    patchIDList_ = -1;
-
-    forAll(nearest, sampleI)
-    {
-        processor_[sampleI] = nearest[sampleI].second().second();
-
-        if (nearest[sampleI].second().second() == Pstream::myProcNo())
-        {
-            // Store the face to sample
-            faceList_[sampleI] = nearest[sampleI].first().index();
-            const label facei = faceList_[sampleI];
-            if (facei != -1)
-            {
-                processor_[sampleI] = Pstream::myProcNo();
-                patchIDList_[sampleI] = bm.whichPatch(facei);
-            }
-        }
-        reduce(processor_[sampleI], maxOp<label>());
-        reduce(patchIDList_[sampleI], maxOp<label>());
-    }
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::patchProbes::patchProbes
@@ -232,7 +53,14 @@ Foam::patchProbes::patchProbes
     const bool readFields
 )
 :
-    probes(name, runTime, dict, loadFromFiles, false)
+    Base
+    (
+        name,
+        runTime,
+        dict,
+        loadFromFiles,
+        readFields
+    )
 {
     if (readFields)
     {
@@ -243,52 +71,13 @@ Foam::patchProbes::patchProbes
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-bool Foam::patchProbes::performAction(unsigned request)
-{
-    if (!pointField::empty() && request && prepare(request))
-    {
-        performAction(scalarFields_, request);
-        performAction(vectorFields_, request);
-        performAction(sphericalTensorFields_, request);
-        performAction(symmTensorFields_, request);
-        performAction(tensorFields_, request);
-
-        performAction(surfaceScalarFields_, request);
-        performAction(surfaceVectorFields_, request);
-        performAction(surfaceSphericalTensorFields_, request);
-        performAction(surfaceSymmTensorFields_, request);
-        performAction(surfaceTensorFields_, request);
-    }
-    return true;
-}
-
-
-bool Foam::patchProbes::execute()
-{
-    if (onExecute_)
-    {
-        return performAction(ACTION_ALL & ~ACTION_WRITE);
-    }
-
-    return true;
-}
-
-
-bool Foam::patchProbes::write()
-{
-    return performAction(ACTION_ALL);
-}
-
-
 bool Foam::patchProbes::read(const dictionary& dict)
 {
-    if (!dict.readIfPresent("patches", patchNames_))
+    if (!(Probes::read(dict) && prober_.read(dict)))
     {
-        patchNames_.resize(1);
-        patchNames_.first() = dict.get<word>("patch");
+        return false;
     }
-
-    return probes::read(dict);
+    return true;
 }
 
 

src/sampling/probes/patchProbes.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2022 OpenCFD Ltd.
+    Copyright (C) 2016-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,52 +28,71 @@ Class
     Foam::patchProbes
 
 Description
-    Set of locations to sample at patches
+    Probes the specified points on specified patches. The points get snapped
+    onto the nearest point on the nearest face of the specified patch, and the
+    sampling is actioned on the snapped locations.
 
-    Call write() to sample and write files.
-    - find nearest location on nearest face
-    - update *this with location (so header contains 'snapped' locations
-    - use *this as the sampling location
-
-    Example of function object specification:
+Usage
+    Minimal example by using \c system/controlDict.functions:
     \verbatim
     patchProbes
     {
+        // Mandatory entries
         type            patchProbes;
         libs            (sampling);
 
-        // Name of the directory for probe data
-        name            patchProbes;
+        fields          (<wordRes>);
+        probeLocations  (<vectorList>);
+        patches         (<wordRes>);    // or patch  <word>;
 
-        // Patches to sample (wildcards allowed)
-        patches         (".*inl.*");
+        // Optional entries
+        verbose         <bool>;
+        sampleOnExecute <bool>;
+        fixedLocations  <bool>;
+        includeOutOfBounds <bool>;
+        interpolationScheme <word>;
 
-        // Write at same frequency as fields
-        writeControl    writeTime;
-        writeInterval   1;
-
-        // Fields to be probed
-        fields          (p U);
-
-        // Locations to probe. These get snapped onto the nearest point
-        // on the selected patches
-        probeLocations
-        (
-            ( -100 0 0.01 )      // at inlet
-        );
+        ...
     }
     \endverbatim
 
+    where the entries mean:
+    \table
+      Property     | Description                        | Type | Reqd | Deflt
+      type         | Type name: patchProbes             | word | yes  | -
+      libs         | Library name: sampling             | word | yes  | -
+      fields       | Names of the fields to be probed   | wordRes | yes | -
+      probeLocations | Locations of the probes          | vectorField | yes | -
+      patches      | Patches to sample (wildcards allowed) | wordRes | yes | -
+      verbose      | Enable/disable verbose output         | bool | no | false
+      sampleOnExecute | Sample on execution and store results | bool | no <!--
+                   -->                                               | false
+      fixedLocations | Do not recalculate cells if mesh moves | bool | no | true
+      includeOutOfBounds | Include out-of-bounds locations | bool | no | true
+      interpolationScheme | Scheme to obtain values at the points | word <!--
+                   -->                                            | no | cell
+    \endtable
+
+    The inherited entries are elaborated in:
+     - \link fvMeshFunctionObject.H \endlink
+     - \link Probes.H \endlink
+     - \link patchProber.H \endlink
+     - \link prober.H \endlink
+
+Note
+  - The \c includeOutOfBounds filters out points that haven't been found.
+    Default is to include them (with value \c -VGREAT).
+
 SourceFiles
     patchProbes.C
-    patchProbesTemplates.C
 
 \*---------------------------------------------------------------------------*/
 
 #ifndef Foam_patchProbes_H
 #define Foam_patchProbes_H
 
-#include "probes.H"
+#include "Probes.H"
+#include "patchProber.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -81,57 +100,17 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                          Class patchProbes Declaration
+                        Class patchProbes Declaration
 \*---------------------------------------------------------------------------*/
 
 class patchProbes
 :
-    public probes
+    public Probes<patchProber>
 {
-protected:
+    // Private Data
 
-    // Protected Data
-
-        //- Patches to sample
-        wordRes patchNames_;
-
-
-    // Protected Member Functions
-
-        //- Find elements containing patchProbes
-        virtual void findElements(const fvMesh& mesh);  // override
-
-
-private:
-
-    // Private Member Functions
-
-        //- Write field values
-        template<class Type>
-        void writeValues
-        (
-            const word& fieldName,
-            const Field<Type>& values,
-            const scalar timeValue
-        );
-
-        //- Sample and store/write applicable volume/surface fields
-        template<class GeoField>
-        void performAction
-        (
-            const fieldGroup<GeoField>& fieldNames,  /* must be sorted */
-            unsigned request
-        );
-
-        //- Perform sampling action with store/write
-        bool performAction(unsigned request);
-
-
-        //- No copy construct
-        patchProbes(const patchProbes&) = delete;
-
-        //- No copy assignment
-        void operator=(const patchProbes&) = delete;
+        //- Use simpler synonym for the base type
+        using Base = Probes<patchProber>;
 
 
 public:
@@ -142,7 +121,7 @@ public:
 
     // Constructors
 
-        //- Construct from Time and dictionary
+        //- Construct from name, Time and dictionary
         patchProbes
         (
             const word& name,
@@ -152,54 +131,26 @@ public:
             const bool readFields = true
         );
 
+
     //- Destructor
     virtual ~patchProbes() = default;
 
 
     // Member Functions
 
-        //- Sample and store result if the sampleOnExecute is enabled.
-        virtual bool execute();
+        //- Bring Base::prober into this class's public scope.
+        using Base::prober;
 
-        //- Sample and write
-        virtual bool write();
-
-        //- Read
+        //- Read the settings from the dictionary
         virtual bool read(const dictionary&);
-
-
-    // Sampling
-
-        //- Sample a volume field at all locations
-        template<class Type>
-        tmp<Field<Type>> sample(const VolumeField<Type>&) const;
-
-        //- Sample a surface field at all locations
-        template<class Type>
-        tmp<Field<Type>> sample(const SurfaceField<Type>&) const;
-
-        //- Sample a single field on all sample locations
-        template<class Type>
-        tmp<Field<Type>> sample(const word& fieldName) const;
-
-        //- Sample a surface field at all locations
-        template<class Type>
-        tmp<Field<Type>> sampleSurfaceField(const word& fieldName) const;
 };
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#ifdef NoRepository
-    #include "patchProbesTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //

src/sampling/probes/probers/internalProber/internalProber.C

@@ -0,0 +1,188 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "internalProber.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(internalProber, 0);
+}
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void Foam::internalProber::findElements(const fvMesh& mesh)
+{
+    DebugInfo<< "internalProber: resetting sample locations" << endl;
+
+    const pointField& probeLocations = this->probeLocations();
+
+    cellIds_.resize_nocopy(probeLocations.size());
+    faceIds_.resize_nocopy(probeLocations.size());
+    procIds_.resize_nocopy(probeLocations.size());
+    procIds_ = -1;
+
+    forAll(probeLocations, probei)
+    {
+        const point& location = probeLocations[probei];
+
+        const label celli = mesh.findCell(location);
+
+        cellIds_[probei] = celli;
+
+        if (celli != -1)
+        {
+            const labelList& cellFaces = mesh.cells()[celli];
+            const vector& cellCentre = mesh.cellCentres()[celli];
+            scalar minDistance = GREAT;
+            label minFaceID = -1;
+            forAll(cellFaces, i)
+            {
+                label facei = cellFaces[i];
+                vector dist = mesh.faceCentres()[facei] - cellCentre;
+                if (mag(dist) < minDistance)
+                {
+                    minDistance = mag(dist);
+                    minFaceID = facei;
+                }
+            }
+            faceIds_[probei] = minFaceID;
+        }
+        else
+        {
+            faceIds_[probei] = -1;
+        }
+
+        if (debug && (cellIds_[probei] != -1 || faceIds_[probei] != -1))
+        {
+            Pout<< "internalProber : found point " << location
+                << " in cell " << cellIds_[probei]
+                << " and face " << faceIds_[probei] << endl;
+        }
+    }
+
+
+    // Check if all probes have been found.
+    forAll(cellIds_, probei)
+    {
+        const point& location = probeLocations[probei];
+        label celli = cellIds_[probei];
+        label facei = faceIds_[probei];
+
+        procIds_[probei] = (celli != -1 ? Pstream::myProcNo() : -1);
+
+        // Check at least one processor with cell.
+        reduce(celli, maxOp<label>());
+        reduce(facei, maxOp<label>());
+        reduce(procIds_[probei], maxOp<label>());
+
+        if (celli == -1)
+        {
+            if (Pstream::master())
+            {
+                WarningInFunction
+                    << "Did not find location " << location
+                    << " in any cell. Skipping location." << endl;
+            }
+        }
+        else if (facei == -1)
+        {
+            if (Pstream::master())
+            {
+                WarningInFunction
+                    << "Did not find location " << location
+                    << " in any face. Skipping location." << endl;
+            }
+        }
+        else
+        {
+            // Make sure location not on two domains.
+            if (cellIds_[probei] != -1 && cellIds_[probei] != celli)
+            {
+                WarningInFunction
+                    << "Location " << location
+                    << " seems to be on multiple domains:"
+                    << " cell " << cellIds_[probei]
+                    << " on my domain " << Pstream::myProcNo()
+                    << " and cell " << celli << " on some other domain."
+                    << nl
+                    << "This might happen if the probe location is on"
+                    << " a processor patch. Change the location slightly"
+                    << " to prevent this." << endl;
+            }
+
+            if (faceIds_[probei] != -1 && faceIds_[probei] != facei)
+            {
+                WarningInFunction
+                    << "Location " << location
+                    << " seems to be on multiple domains:"
+                    << " cell " << faceIds_[probei]
+                    << " on my domain " << Pstream::myProcNo()
+                    << " and face " << facei << " on some other domain."
+                    << nl
+                    << "This might happen if the probe location is on"
+                    << " a processor patch. Change the location slightly"
+                    << " to prevent this." << endl;
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::internalProber::internalProber
+(
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    prober(mesh, dict)
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::internalProber::read(const dictionary& dict)
+{
+    if (!prober::read(dict))
+    {
+        return false;
+    }
+
+    // Initialise cells to sample from supplied locations
+    findElements(mesh_);
+
+    return true;
+}
+
+
+// ************************************************************************* //

src/sampling/probes/probers/internalProber/internalProber.H

@@ -0,0 +1,140 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::internalProber
+
+Description
+    A utility class for probing field values at specified point locations
+    within an \c fvMesh.
+
+    The \c internalProber stores a list of 3D point coordinates and
+    determines the corresponding mesh elements (cells or faces) that contain
+    these points. It provides methods to sample volume or surface fields at
+    the stored locations, with support for fixed or mesh-moving point
+    coordinates.
+
+    Features include:
+      - Reading probe locations and settings from a dictionary
+      - Support for fixed or moving locations (for dynamic mesh cases)
+      - Optional inclusion of points that lie outside of the mesh domain
+      - Selection of interpolation/sampling schemes for fixed locations
+      - Sampling of volume and surface fields by name or by direct reference
+      - Automatic update of element mapping when the mesh changes or moves
+
+SourceFiles
+    internalProber.C
+    internalProberTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_internalProber_H
+#define Foam_internalProber_H
+
+#include "prober.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class internalProber Declaration
+\*---------------------------------------------------------------------------*/
+
+class internalProber
+:
+    public prober
+{
+protected:
+
+    // Protected Member Functions
+
+        //- Find cells and faces containing probes
+        virtual void findElements(const fvMesh& mesh);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("internalProber");
+
+
+    // Constructors
+
+        //- Construct from Time and dictionary
+        internalProber
+        (
+            const fvMesh& mesh,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~internalProber() = default;
+
+
+    // Member Functions
+
+    // Sampling
+
+        //- Sample a volume field at all locations
+        template<class Type>
+        tmp<Field<Type>> sample(const VolumeField<Type>&) const;
+
+        //- Sample a surface field at all locations
+        template<class Type>
+        tmp<Field<Type>> sample(const SurfaceField<Type>&) const;
+
+        //- Sample a volume field at all locations
+        template<class Type>
+        tmp<Field<Type>> sample(const word& fieldName) const;
+
+        //- Sample a surface field at all locations
+        template<class Type>
+        tmp<Field<Type>> sampleSurfaceField(const word& fieldName) const;
+
+
+    // I-O
+
+        //- Read the settings dictionary
+        virtual bool read(const dictionary&);
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "internalProberTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/sampling/probes/probers/internalProber/internalProberTemplates.C

@@ -0,0 +1,122 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "internalProber.H"
+#include "interpolation.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::internalProber::sample(const VolumeField<Type>& vField) const
+{
+    const Type unsetVal(-VGREAT*pTraits<Type>::one);
+
+    auto tvalues = tmp<Field<Type>>::New(Field<Type>(this->size(), unsetVal));
+    auto& values = tvalues.ref();
+
+    if (fixedLocations_)
+    {
+        autoPtr<interpolation<Type>> interpPtr
+        (
+            interpolation<Type>::New(samplePointScheme_, vField)
+        );
+
+        const pointField& probeLocations = this->probeLocations();
+        forAll(probeLocations, probei)
+        {
+            if (cellIds_[probei] >= 0)
+            {
+                const vector& position = probeLocations[probei];
+
+                values[probei] = interpPtr().interpolate
+                (
+                    position,
+                    cellIds_[probei],
+                    -1
+                );
+            }
+        }
+    }
+    else
+    {
+        forAll(*this, probei)
+        {
+            if (cellIds_[probei] >= 0)
+            {
+                values[probei] = vField[cellIds_[probei]];
+            }
+        }
+    }
+
+    Pstream::listCombineReduce(values, isNotEqOp<Type>());
+
+    return tvalues;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::internalProber::sample(const SurfaceField<Type>& sField) const
+{
+    const Type unsetVal(-VGREAT*pTraits<Type>::one);
+
+    auto tvalues = tmp<Field<Type>>::New(Field<Type>(this->size(), unsetVal));
+    auto& values = tvalues.ref();
+
+    const pointField& probeLocations = this->probeLocations();
+    forAll(probeLocations, probei)
+    {
+        if (faceIds_[probei] >= 0)
+        {
+            values[probei] = sField[faceIds_[probei]];
+        }
+    }
+
+    Pstream::listCombineReduce(values, isNotEqOp<Type>());
+
+    return tvalues;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::internalProber::sample(const word& fieldName) const
+{
+    return sample(mesh_.lookupObject<VolumeField<Type>>(fieldName));
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::internalProber::sampleSurfaceField(const word& fieldName) const
+{
+    return sample(mesh_.lookupObject<SurfaceField<Type>>(fieldName));
+}
+
+
+// ************************************************************************* //

src/sampling/probes/probers/patchProber/patchProber.C

@@ -0,0 +1,252 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2016-2022 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 "patchProber.H"
+#include "mappedPatchBase.H"
+#include "treeBoundBox.H"
+#include "treeDataFace.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(patchProber, 0);
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::patchProber::findElements(const fvMesh& mesh)
+{
+    (void)mesh.tetBasePtIs();
+
+    const polyBoundaryMesh& bm = mesh.boundaryMesh();
+
+    // All the info for nearest. Construct to miss
+    List<mappedPatchBase::nearInfo> nearest(this->size());
+
+    patchIDs_ = bm.patchSet(patchNames_).sortedToc();
+
+    label nFaces = 0;
+    forAll(patchIDs_, i)
+    {
+        nFaces += bm[patchIDs_[i]].size();
+    }
+
+    if (nFaces > 0)
+    {
+        // Collect mesh faces and bounding box
+        labelList bndFaces(nFaces);
+        treeBoundBox overallBb;
+
+        nFaces = 0;
+        forAll(patchIDs_, i)
+        {
+            const polyPatch& pp = bm[patchIDs_[i]];
+            forAll(pp, i)
+            {
+                bndFaces[nFaces++] = pp.start()+i;
+                const face& f = pp[i];
+
+                // Without reduction.
+                overallBb.add(pp.points(), f);
+            }
+        }
+
+        Random rndGen(123456);
+        overallBb.inflate(rndGen, 1e-4, ROOTVSMALL);
+
+
+        const indexedOctree<treeDataFace> boundaryTree
+        (
+            treeDataFace(mesh, bndFaces),  // patch faces only
+
+            overallBb,      // overall search domain
+            8,              // maxLevel
+            10,             // leafsize
+            3.0             // duplicity
+        );
+
+        forAll(probeLocations(), probei)
+        {
+            const auto& treeData = boundaryTree.shapes();
+            const point sample = probeLocations()[probei];
+
+            pointIndexHit info = boundaryTree.findNearest
+            (
+                sample,
+                Foam::sqr(boundaryTree.bb().mag())
+            );
+
+            if (!info.hit())
+            {
+                info = boundaryTree.findNearest(sample, Foam::sqr(GREAT));
+            }
+
+            const label facei = treeData.objectIndex(info.index());
+
+            const label patchi = bm.whichPatch(facei);
+
+            if (isA<emptyPolyPatch>(bm[patchi]))
+            {
+                WarningInFunction
+                    << " The sample point: " << sample
+                    << " belongs to " << patchi
+                    << " which is an empty patch. This is not permitted. "
+                    << " This sample will not be included "
+                    << endl;
+            }
+            else if (info.hit())
+            {
+                // Note: do we store the face centre or the actual nearest?
+                // We interpolate using the faceI only though (no
+                // interpolation) so it does not actually matter much, just for
+                // the location written to the header.
+
+                //const point& facePt = mesh.faceCentres()[faceI];
+                const point& facePt = info.point();
+
+                mappedPatchBase::nearInfo sampleInfo;
+
+                sampleInfo.first() = pointIndexHit(true, facePt, facei);
+
+                sampleInfo.second().first() = facePt.distSqr(sample);
+                sampleInfo.second().second() = Pstream::myProcNo();
+
+                nearest[probei] = sampleInfo;
+            }
+        }
+    }
+
+
+    // Find nearest - globally consistent
+    Pstream::listCombineReduce(nearest, mappedPatchBase::nearestEqOp());
+
+    oldPoints_.resize(this->size());
+
+    pointField& probeLocations = this->probeLocations();
+    // Update actual probe locations and store old ones
+    forAll(nearest, samplei)
+    {
+        oldPoints_[samplei] = probeLocations[samplei];
+        probeLocations[samplei] = nearest[samplei].first().point();
+    }
+
+    if (debug)
+    {
+        InfoInFunction << nl;
+        forAll(nearest, samplei)
+        {
+            label proci = nearest[samplei].second().second();
+            label locali = nearest[samplei].first().index();
+
+            Info<< "    " << samplei << " coord:"<< probeLocations[samplei]
+                << " found on processor:" << proci
+                << " in local face:" << locali
+                << " with location:" << nearest[samplei].first().point()
+                << endl;
+        }
+    }
+
+    // Extract any local faces to sample:
+    // - operator[] : actual point to sample (=nearest point on patch)
+    // - oldPoints_ : original provided point (might be anywhere in the mesh)
+    // - cellIds_   : cells, not used
+    // - faceIds_      : faces (now patch faces)
+    // - patchIds_   : patch corresponding to faceList
+    // - procIds_     : processor
+    cellIds_.resize_nocopy(nearest.size());
+    cellIds_ = -1;
+
+    faceIds_.resize_nocopy(nearest.size());
+    faceIds_ = -1;
+
+    procIds_.resize_nocopy(nearest.size());
+    procIds_ = -1;
+
+    patchIds_.resize_nocopy(nearest.size());
+    patchIds_ = -1;
+
+    forAll(nearest, sampleI)
+    {
+        procIds_[sampleI] = nearest[sampleI].second().second();
+
+        if (nearest[sampleI].second().second() == Pstream::myProcNo())
+        {
+            // Store the face to sample
+            faceIds_[sampleI] = nearest[sampleI].first().index();
+            const label facei = faceIds_[sampleI];
+            if (facei != -1)
+            {
+                procIds_[sampleI] = Pstream::myProcNo();
+                patchIds_[sampleI] = bm.whichPatch(facei);
+            }
+        }
+        reduce(procIds_[sampleI], maxOp<label>());
+        reduce(patchIds_[sampleI], maxOp<label>());
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::patchProber::patchProber
+(
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    prober(mesh, dict)
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::patchProber::read(const dictionary& dict)
+{
+    if (!prober::read(dict))
+    {
+        return false;
+    }
+
+    if (!dict.readIfPresent("patches", patchNames_))
+    {
+        patchNames_.resize(1);
+        patchNames_.first() = dict.get<word>("patch");
+    }
+
+    // Initialise cells to sample from supplied locations
+    findElements(mesh_);
+
+    return true;
+}
+
+
+// ************************************************************************* //

src/sampling/probes/probers/patchProber/patchProber.H

@@ -0,0 +1,152 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::patchProber
+
+Description
+    Utility class for probing specified points on user-selected boundary
+    patches. The input points are projected onto the nearest point of the
+    nearest face on the specified patch, ensuring sampling occurs at valid
+    patch locations.
+
+    The patchProber enables sampling of both volume and surface fields
+    at these snapped locations. Patch selection is controlled via patch names or
+    indices, and the class provides runtime selection and dictionary-driven
+    configuration.
+
+    Typical usage involves specifying patch names and probe locations in a
+    dictionary, after which the class manages the mapping and sampling
+    operations.
+
+SourceFiles
+    patchProber.C
+    patchProberTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_patchProber_H
+#define Foam_patchProber_H
+
+#include "prober.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class patchProber Declaration
+\*---------------------------------------------------------------------------*/
+
+class patchProber
+:
+    public prober
+{
+protected:
+
+    // Protected Data
+
+        //- Names of the patches to sample
+        wordRes patchNames_;
+
+        //- Index of the patches to sample
+        labelList patchIDs_;
+
+
+    // Protected Member Functions
+
+        //- Find elements containing patchProber
+        virtual void findElements(const fvMesh& mesh);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("patchProber");
+
+
+    // Constructors
+
+        //- Construct from Time and dictionary
+        patchProber
+        (
+            const fvMesh& mesh,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~patchProber() = default;
+
+
+    // Member Functions
+
+    // Access
+
+        //- Return the index of the patches to sample
+        const labelList& patchIDs() const noexcept { return patchIDs_; }
+
+
+    // Sampling
+
+        //- Sample a volume field at all locations
+        template<class Type>
+        tmp<Field<Type>> sample(const VolumeField<Type>&) const;
+
+        //- Sample a surface field at all locations
+        template<class Type>
+        tmp<Field<Type>> sample(const SurfaceField<Type>&) const;
+
+        //- Sample a volume field at all locations
+        template<class Type>
+        tmp<Field<Type>> sample(const word& fieldName) const;
+
+        //- Sample a surface field at all locations
+        template<class Type>
+        tmp<Field<Type>> sampleSurfaceField(const word& fieldName) const;
+
+
+    // I-O
+
+        //- Read the settings dictionary
+        virtual bool read(const dictionary&);
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "patchProberTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/sampling/probes/probers/patchProber/patchProberTemplates.C

@@ -5,8 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2021-2022 OpenCFD Ltd.
+    Copyright (C) 2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -26,70 +25,15 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "patchProbes.H"
+#include "patchProber.H"
 #include "volFields.H"
 #include "surfaceFields.H"
-#include "IOmanip.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-template<class Type>
-void Foam::patchProbes::writeValues
-(
-    const word& fieldName,
-    const Field<Type>& values,
-    const scalar timeValue
-)
-{
-    if (Pstream::master())
-    {
-        const unsigned int w = IOstream::defaultPrecision() + 7;
-        OFstream& os = *probeFilePtrs_[fieldName];
-
-        os  << setw(w) << timeValue;
-
-        for (const auto& v : values)
-        {
-            os  << ' ' << setw(w) << v;
-        }
-        os  << endl;
-    }
-}
-
-
-template<class GeoField>
-void Foam::patchProbes::performAction
-(
-    const fieldGroup<GeoField>& fieldNames,
-    unsigned request
-)
-{
-    for (const word& fieldName : fieldNames)
-    {
-        tmp<GeoField> tfield = getOrLoadField<GeoField>(fieldName);
-
-        if (tfield)
-        {
-            const auto& field = tfield();
-            const scalar timeValue = field.time().timeOutputValue();
-
-            Field<typename GeoField::value_type> values(sample(field));
-
-            this->storeResults(fieldName, values);
-            if (request & ACTION_WRITE)
-            {
-                this->writeValues(fieldName, values, timeValue);
-            }
-        }
-    }
-}
-
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
 Foam::tmp<Foam::Field<Type>>
-Foam::patchProbes::sample(const VolumeField<Type>& vField) const
+Foam::patchProber::sample(const VolumeField<Type>& vField) const
 {
     const Type unsetVal(-VGREAT*pTraits<Type>::one);
 
@@ -101,7 +45,7 @@ Foam::patchProbes::sample(const VolumeField<Type>& vField) const
 
     forAll(*this, probei)
     {
-        label facei = faceList_[probei];
+        label facei = faceIds_[probei];
 
         if (facei >= 0)
         {
@@ -119,7 +63,7 @@ Foam::patchProbes::sample(const VolumeField<Type>& vField) const
 
 template<class Type>
 Foam::tmp<Foam::Field<Type>>
-Foam::patchProbes::sample(const SurfaceField<Type>& sField) const
+Foam::patchProber::sample(const SurfaceField<Type>& sField) const
 {
     const Type unsetVal(-VGREAT*pTraits<Type>::one);
 
@@ -131,7 +75,7 @@ Foam::patchProbes::sample(const SurfaceField<Type>& sField) const
 
     forAll(*this, probei)
     {
-        label facei = faceList_[probei];
+        label facei = faceIds_[probei];
 
         if (facei >= 0)
         {
@@ -149,7 +93,7 @@ Foam::patchProbes::sample(const SurfaceField<Type>& sField) const
 
 template<class Type>
 Foam::tmp<Foam::Field<Type>>
-Foam::patchProbes::sample(const word& fieldName) const
+Foam::patchProber::sample(const word& fieldName) const
 {
     return sample(mesh_.lookupObject<VolumeField<Type>>(fieldName));
 }
@@ -157,7 +101,7 @@ Foam::patchProbes::sample(const word& fieldName) const
 
 template<class Type>
 Foam::tmp<Foam::Field<Type>>
-Foam::patchProbes::sampleSurfaceField(const word& fieldName) const
+Foam::patchProber::sampleSurfaceField(const word& fieldName) const
 {
     return sample(mesh_.lookupObject<SurfaceField<Type>>(fieldName));
 }

src/sampling/probes/probers/prober.C

@@ -0,0 +1,188 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "prober.H"
+#include "mapPolyMesh.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(prober, 0);
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::prober::prober
+(
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    mesh_(mesh),
+    samplePointScheme_("cell")
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::prober::read(const dictionary& dict)
+{
+    dict.readEntry("probeLocations", probes_);
+
+    if (probes_.empty())
+    {
+        FatalIOErrorInFunction(dict)
+            << "Empty 'probeLocations' list."
+            << exit(FatalIOError);
+    }
+
+    fixedLocations_ = dict.getOrDefault<bool>("fixedLocations", true);
+    includeOutOfBounds_ = dict.getOrDefault<bool>("includeOutOfBounds", true);
+
+    if (dict.readIfPresent("interpolationScheme", samplePointScheme_))
+    {
+        if (!fixedLocations_ && samplePointScheme_ != "cell")
+        {
+            WarningInFunction
+                << "Only cell interpolation can be applied when "
+                << "not using fixedLocations. InterpolationScheme "
+                << "entry will be ignored"
+                << endl;
+        }
+    }
+
+    return true;
+}
+
+
+void Foam::prober::updateMesh(const mapPolyMesh& mpm)
+{
+    DebugInfo<< "probes: updateMesh" << endl;
+
+    if (&mpm.mesh() != &mesh_)
+    {
+        return;
+    }
+
+    if (fixedLocations_)
+    {
+        this->findElements(mesh_);
+    }
+    else
+    {
+        DebugInfo<< "probes: remapping sample locations" << endl;
+
+        // 1. Update cells
+        {
+            DynamicList<label> elems(cellIds_.size());
+
+            const labelList& reverseMap = mpm.reverseCellMap();
+            forAll(cellIds_, i)
+            {
+                label celli = cellIds_[i];
+                if (celli != -1)
+                {
+                    label newCelli = reverseMap[celli];
+                    if (newCelli == -1)
+                    {
+                        // cell removed
+                    }
+                    else if (newCelli < -1)
+                    {
+                        // cell merged
+                        elems.append(-newCelli - 2);
+                    }
+                    else
+                    {
+                        // valid new cell
+                        elems.append(newCelli);
+                    }
+                }
+                else
+                {
+                    // Keep -1 elements so the size stays the same
+                    elems.append(-1);
+                }
+            }
+
+            cellIds_.transfer(elems);
+        }
+
+        // 2. Update faces
+        {
+            DynamicList<label> elems(faceIds_.size());
+
+            const labelList& reverseMap = mpm.reverseFaceMap();
+            for (const label facei : faceIds_)
+            {
+                if (facei != -1)
+                {
+                    label newFacei = reverseMap[facei];
+                    if (newFacei == -1)
+                    {
+                        // face removed
+                    }
+                    else if (newFacei < -1)
+                    {
+                        // face merged
+                        elems.append(-newFacei - 2);
+                    }
+                    else
+                    {
+                        // valid new face
+                        elems.append(newFacei);
+                    }
+                }
+                else
+                {
+                    // Keep -1 elements
+                    elems.append(-1);
+                }
+            }
+
+            faceIds_.transfer(elems);
+        }
+    }
+}
+
+
+void Foam::prober::movePoints(const polyMesh& mesh)
+{
+    DebugInfo<< "probes: movePoints" << endl;
+
+    if (fixedLocations_ && &mesh == &mesh_)
+    {
+        this->findElements(mesh_);
+    }
+}
+
+
+// ************************************************************************* //

src/sampling/probes/probers/prober.H

@@ -0,0 +1,220 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::patchProber
+
+Description
+    Base class for sampling fields at specified internal and boundary locations.
+
+SourceFiles
+    prober.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_prober_H
+#define Foam_prober_H
+
+#include "fvMesh.H"
+#include "pointField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class prober Declaration
+\*---------------------------------------------------------------------------*/
+
+class prober
+{
+protected:
+
+    template<class T>
+    struct isNotEqOp
+    {
+        void operator()(T& x, const T& y) const
+        {
+            const T unsetVal(-VGREAT*pTraits<T>::one);
+
+            if (x != unsetVal)
+            {
+                // Keep x.
+
+                // Note: should check for y != unsetVal but multiple sample cells
+                // already handled in read().
+            }
+            else
+            {
+                // x is not set. y might be.
+                x = y;
+            }
+        }
+    };
+
+
+    // Protected Data
+
+        //- Const reference to the mesh
+        const fvMesh& mesh_;
+
+        //- Fixed locations (default: true)
+        //  Note: set to false for moving mesh calculations where locations
+        //        should move with the mesh
+        bool fixedLocations_;
+
+        //- Include probes that were not found (default: true)
+        bool includeOutOfBounds_;
+
+        //- Interpolation/sample scheme to obtain values at the points
+        //  Note: only possible when fixedLocations_ is true
+        word samplePointScheme_;
+
+
+    // Calculated
+
+        //- Probe locations
+        pointField probes_;
+
+        //- Cells to be probed (obtained from the locations)
+        labelList cellIds_;
+
+        //- Faces to be probed
+        labelList faceIds_;
+
+        //- Processor holding the cell or face (-1 if point not found
+        //- on any processor)
+        labelList procIds_;
+
+        //- Patch IDs on which the new probes are located
+        labelList patchIds_;
+
+        //- Original probes location
+        pointField oldPoints_;
+
+
+    // Protected Member Functions
+
+        //- Find cells and faces containing probes
+        virtual void findElements(const fvMesh& mesh) = 0;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("prober");
+
+
+    // Generated Methods
+
+        //- No copy construct
+        prober(const prober&) = delete;
+
+        //- No copy assignment
+        void operator=(const prober&) = delete;
+
+
+    // Constructors
+
+        //- Construct from Time and dictionary
+        prober
+        (
+            const fvMesh& mesh,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~prober() = default;
+
+
+    // Member Functions
+
+    // Access
+
+        //- Return true if no probe locations
+        bool empty() const { return probes_.empty(); }
+
+        //- Return number of probe locations
+        label size() const { return probes_.size(); }
+
+        //- Return true if fixed locations
+        bool fixedLocations() const { return fixedLocations_; }
+
+        //- Return true if include out of bounds probes
+        bool includeOutOfBounds() const { return includeOutOfBounds_; }
+
+        //- Return the interpolation scheme to obtain values at the points
+        //  Note: only possible when fixedLocations_ is true
+        const word& samplePointScheme() const { return samplePointScheme_; }
+
+        //- Return const reference to the probe locations
+        const pointField& probeLocations() const { return probes_; }
+
+        //- Return reference to the probe locations
+        pointField& probeLocations() { return probes_; }
+
+        //- Return the location of probe i
+        const point& probe(const label i) const { return probes_[i]; }
+
+        //- Cells to be probed (obtained from the locations)
+        const labelList& elements() const { return cellIds_; }
+
+        //- Return const reference to the faces to be probed
+        const labelList& faces() const { return faceIds_; }
+
+        //- Return const reference to the processor list
+        const labelList& processors() const { return procIds_; }
+
+        //- Return const reference to the patch ID list
+        const labelList& patchIDList() const noexcept { return patchIds_; }
+
+        //- Return const reference to the original probe locations
+        const pointField& oldPoints() const noexcept { return oldPoints_; }
+
+
+    // I-O
+
+        //- Read the settings dictionary
+        virtual bool read(const dictionary&);
+
+        //- Update for changes of mesh
+        virtual void updateMesh(const mapPolyMesh&);
+
+        //- Update for changes of mesh
+        virtual void movePoints(const polyMesh&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/sampling/probes/probes.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2015-2023 OpenCFD Ltd.
+    Copyright (C) 2015-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,12 +27,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "probes.H"
-#include "dictionary.H"
-#include "volFields.H"
-#include "surfaceFields.H"
-#include "Time.H"
-#include "IOmanip.H"
-#include "mapPolyMesh.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -40,7 +34,6 @@ License
 namespace Foam
 {
     defineTypeNameAndDebug(probes, 0);
-
     addToRunTimeSelectionTable
     (
         functionObject,
@@ -49,315 +42,6 @@ namespace Foam
     );
 }
 
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-void Foam::probes::createProbeFiles(const wordList& fieldNames)
-{
-    // Open new output streams
-
-    bool needsNewFiles = false;
-    for (const word& fieldName : fieldNames)
-    {
-        if (!probeFilePtrs_.found(fieldName))
-        {
-            needsNewFiles = true;
-            break;
-        }
-    }
-
-    if (needsNewFiles && Pstream::master())
-    {
-        DebugInfo
-            << "Probing fields: " << fieldNames << nl
-            << "Probing locations: " << *this << nl
-            << endl;
-
-        // Put in undecomposed case
-        // (Note: gives problems for distributed data running)
-
-        fileName probeDir
-        (
-            mesh_.time().globalPath()
-          / functionObject::outputPrefix
-          / name()/mesh_.regionName()
-            // Use startTime as the instance for output files
-          / mesh_.time().timeName(mesh_.time().startTime().value())
-        );
-        probeDir.clean();  // Remove unneeded ".."
-
-        // Create directory if needed
-        Foam::mkDir(probeDir);
-
-        for (const word& fieldName : fieldNames)
-        {
-            if (probeFilePtrs_.found(fieldName))
-            {
-                // Safety
-                continue;
-            }
-
-            auto osPtr = autoPtr<OFstream>::New(probeDir/fieldName);
-            auto& os = *osPtr;
-
-            probeFilePtrs_.insert(fieldName, osPtr);
-
-            DebugInfo<< "open probe stream: " << os.name() << endl;
-
-            const unsigned int width(IOstream::defaultPrecision() + 7);
-            os.setf(std::ios_base::left);
-
-            forAll(*this, probei)
-            {
-                os  << "# Probe " << probei << ' ' << operator[](probei);
-
-                if (processor_[probei] == -1)
-                {
-                    os  << "  # Not Found";
-                }
-                // Only for patchProbes
-                else if (probei < patchIDList_.size())
-                {
-                    const label patchi = patchIDList_[probei];
-                    if (patchi != -1)
-                    {
-                        const polyBoundaryMesh& bm = mesh_.boundaryMesh();
-                        if
-                        (
-                            patchi < bm.nNonProcessor()
-                         || processor_[probei] == Pstream::myProcNo()
-                        )
-                        {
-                            os  << " at patch " << bm[patchi].name();
-                        }
-                        os  << " with a distance of "
-                            << mag(operator[](probei)-oldPoints_[probei])
-                            << " m to the original point "
-                            << oldPoints_[probei];
-                    }
-                }
-
-                os  << nl;
-            }
-
-            os  << setw(width) << "# Time";
-
-            forAll(*this, probei)
-            {
-                if (includeOutOfBounds_ || processor_[probei] != -1)
-                {
-                    os  << ' ' << setw(width) << probei;
-                }
-            }
-            os  << endl;
-        }
-    }
-}
-
-
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-void Foam::probes::findElements(const fvMesh& mesh)
-{
-    DebugInfo<< "probes: resetting sample locations" << endl;
-
-    elementList_.resize_nocopy(pointField::size());
-    faceList_.resize_nocopy(pointField::size());
-    processor_.resize_nocopy(pointField::size());
-    processor_ = -1;
-
-    forAll(*this, probei)
-    {
-        const point& location = (*this)[probei];
-
-        const label celli = mesh.findCell(location);
-
-        elementList_[probei] = celli;
-
-        if (celli != -1)
-        {
-            const labelList& cellFaces = mesh.cells()[celli];
-            const vector& cellCentre = mesh.cellCentres()[celli];
-            scalar minDistance = GREAT;
-            label minFaceID = -1;
-            forAll(cellFaces, i)
-            {
-                label facei = cellFaces[i];
-                vector dist = mesh.faceCentres()[facei] - cellCentre;
-                if (mag(dist) < minDistance)
-                {
-                    minDistance = mag(dist);
-                    minFaceID = facei;
-                }
-            }
-            faceList_[probei] = minFaceID;
-        }
-        else
-        {
-            faceList_[probei] = -1;
-        }
-
-        if (debug && (elementList_[probei] != -1 || faceList_[probei] != -1))
-        {
-            Pout<< "probes : found point " << location
-                << " in cell " << elementList_[probei]
-                << " and face " << faceList_[probei] << endl;
-        }
-    }
-
-
-    // Check if all probes have been found.
-    forAll(elementList_, probei)
-    {
-        const point& location = operator[](probei);
-        label celli = elementList_[probei];
-        label facei = faceList_[probei];
-
-        processor_[probei] = (celli != -1 ? Pstream::myProcNo() : -1);
-
-        // Check at least one processor with cell.
-        reduce(celli, maxOp<label>());
-        reduce(facei, maxOp<label>());
-        reduce(processor_[probei], maxOp<label>());
-
-        if (celli == -1)
-        {
-            if (Pstream::master())
-            {
-                WarningInFunction
-                    << "Did not find location " << location
-                    << " in any cell. Skipping location." << endl;
-            }
-        }
-        else if (facei == -1)
-        {
-            if (Pstream::master())
-            {
-                WarningInFunction
-                    << "Did not find location " << location
-                    << " in any face. Skipping location." << endl;
-            }
-        }
-        else
-        {
-            // Make sure location not on two domains.
-            if (elementList_[probei] != -1 && elementList_[probei] != celli)
-            {
-                WarningInFunction
-                    << "Location " << location
-                    << " seems to be on multiple domains:"
-                    << " cell " << elementList_[probei]
-                    << " on my domain " << Pstream::myProcNo()
-                    << " and cell " << celli << " on some other domain."
-                    << nl
-                    << "This might happen if the probe location is on"
-                    << " a processor patch. Change the location slightly"
-                    << " to prevent this." << endl;
-            }
-
-            if (faceList_[probei] != -1 && faceList_[probei] != facei)
-            {
-                WarningInFunction
-                    << "Location " << location
-                    << " seems to be on multiple domains:"
-                    << " cell " << faceList_[probei]
-                    << " on my domain " << Pstream::myProcNo()
-                    << " and face " << facei << " on some other domain."
-                    << nl
-                    << "This might happen if the probe location is on"
-                    << " a processor patch. Change the location slightly"
-                    << " to prevent this." << endl;
-            }
-        }
-    }
-}
-
-
-Foam::label Foam::probes::prepare(unsigned request)
-{
-    // Prefilter on selection
-    HashTable<wordHashSet> selected =
-    (
-        loadFromFiles_
-      ? IOobjectList(mesh_, mesh_.time().timeName()).classes(fieldSelection_)
-      : mesh_.classes(fieldSelection_)
-    );
-
-    // Classify and count fields
-    label nFields = 0;
-    do
-    {
-        #undef  doLocalCode
-        #define doLocalCode(InputType, Target)                                \
-        {                                                                     \
-            Target.clear();  /* Remove old values */                          \
-            const auto iter = selected.cfind(InputType::typeName);            \
-            if (iter.good())                                                  \
-            {                                                                 \
-                /* Add new (current) values */                                \
-                Target.append(iter.val().sortedToc());                        \
-                nFields += Target.size();                                     \
-            }                                                                 \
-        }
-
-        doLocalCode(volScalarField, scalarFields_);
-        doLocalCode(volVectorField, vectorFields_)
-        doLocalCode(volSphericalTensorField, sphericalTensorFields_);
-        doLocalCode(volSymmTensorField, symmTensorFields_);
-        doLocalCode(volTensorField, tensorFields_);
-
-        doLocalCode(surfaceScalarField, surfaceScalarFields_);
-        doLocalCode(surfaceVectorField, surfaceVectorFields_);
-        doLocalCode(surfaceSphericalTensorField, surfaceSphericalTensorFields_);
-        doLocalCode(surfaceSymmTensorField, surfaceSymmTensorFields_);
-        doLocalCode(surfaceTensorField, surfaceTensorFields_);
-        #undef doLocalCode
-    }
-    while (false);
-
-
-    // Adjust file streams
-    if (Pstream::master())
-    {
-        wordHashSet currentFields(2*nFields);
-        currentFields.insert(scalarFields_);
-        currentFields.insert(vectorFields_);
-        currentFields.insert(sphericalTensorFields_);
-        currentFields.insert(symmTensorFields_);
-        currentFields.insert(tensorFields_);
-
-        currentFields.insert(surfaceScalarFields_);
-        currentFields.insert(surfaceVectorFields_);
-        currentFields.insert(surfaceSphericalTensorFields_);
-        currentFields.insert(surfaceSymmTensorFields_);
-        currentFields.insert(surfaceTensorFields_);
-
-        DebugInfo
-            << "Probing fields: " << currentFields << nl
-            << "Probing locations: " << *this << nl
-            << endl;
-
-        // Close streams for fields that no longer exist
-        forAllIters(probeFilePtrs_, iter)
-        {
-            if (!currentFields.erase(iter.key()))
-            {
-                DebugInfo<< "close probe stream: " << iter()->name() << endl;
-
-                probeFilePtrs_.remove(iter);
-            }
-        }
-
-        if ((request & ACTION_WRITE) && !currentFields.empty())
-        {
-            createProbeFiles(currentFields.sortedToc());
-        }
-    }
-
-    return nFields;
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::probes::probes
@@ -369,15 +53,14 @@ Foam::probes::probes
     const bool readFields
 )
 :
-    functionObjects::fvMeshFunctionObject(name, runTime, dict),
-    pointField(),
-    loadFromFiles_(loadFromFiles),
-    fixedLocations_(true),
-    includeOutOfBounds_(true),
-    verbose_(false),
-    onExecute_(false),
-    fieldSelection_(),
-    samplePointScheme_("cell")
+    Base
+    (
+        name,
+        runTime,
+        dict,
+        loadFromFiles,
+        readFields
+    )
 {
     if (readFields)
     {
@@ -388,184 +71,14 @@ Foam::probes::probes
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-bool Foam::probes::verbose(const bool on) noexcept
-{
-    bool old(verbose_);
-    verbose_ = on;
-    return old;
-}
-
-
 bool Foam::probes::read(const dictionary& dict)
 {
-    dict.readEntry("probeLocations", static_cast<pointField&>(*this));
-    dict.readEntry("fields", fieldSelection_);
-
-    dict.readIfPresent("fixedLocations", fixedLocations_);
-    dict.readIfPresent("includeOutOfBounds", includeOutOfBounds_);
-
-    verbose_ = dict.getOrDefault("verbose", false);
-    onExecute_ = dict.getOrDefault("sampleOnExecute", false);
-
-    if (dict.readIfPresent("interpolationScheme", samplePointScheme_))
+    if (!(Probes::read(dict) && prober_.read(dict)))
     {
-        if (!fixedLocations_ && samplePointScheme_ != "cell")
-        {
-            WarningInFunction
-                << "Only cell interpolation can be applied when "
-                << "not using fixedLocations. InterpolationScheme "
-                << "entry will be ignored"
-                << endl;
-        }
-    }
-
-    // Initialise cells to sample from supplied locations
-    findElements(mesh_);
-
-    // Close old (ununsed) streams
-    prepare(ACTION_NONE);
-
-    return true;
-}
-
-
-bool Foam::probes::performAction(unsigned request)
-{
-    if (!pointField::empty() && request && prepare(request))
-    {
-        performAction(scalarFields_, request);
-        performAction(vectorFields_, request);
-        performAction(sphericalTensorFields_, request);
-        performAction(symmTensorFields_, request);
-        performAction(tensorFields_, request);
-
-        performAction(surfaceScalarFields_, request);
-        performAction(surfaceVectorFields_, request);
-        performAction(surfaceSphericalTensorFields_, request);
-        performAction(surfaceSymmTensorFields_, request);
-        performAction(surfaceTensorFields_, request);
+        return false;
     }
     return true;
 }
 
 
-bool Foam::probes::execute()
-{
-    if (onExecute_)
-    {
-        return performAction(ACTION_ALL & ~ACTION_WRITE);
-    }
-
-    return true;
-}
-
-
-bool Foam::probes::write()
-{
-    return performAction(ACTION_ALL);
-}
-
-
-void Foam::probes::updateMesh(const mapPolyMesh& mpm)
-{
-    DebugInfo<< "probes: updateMesh" << endl;
-
-    if (&mpm.mesh() != &mesh_)
-    {
-        return;
-    }
-
-    if (fixedLocations_)
-    {
-        findElements(mesh_);
-    }
-    else
-    {
-        DebugInfo<< "probes: remapping sample locations" << endl;
-
-        // 1. Update cells
-        {
-            DynamicList<label> elems(elementList_.size());
-
-            const labelList& reverseMap = mpm.reverseCellMap();
-            forAll(elementList_, i)
-            {
-                label celli = elementList_[i];
-                if (celli != -1)
-                {
-                    label newCelli = reverseMap[celli];
-                    if (newCelli == -1)
-                    {
-                        // cell removed
-                    }
-                    else if (newCelli < -1)
-                    {
-                        // cell merged
-                        elems.append(-newCelli - 2);
-                    }
-                    else
-                    {
-                        // valid new cell
-                        elems.append(newCelli);
-                    }
-                }
-                else
-                {
-                    // Keep -1 elements so the size stays the same
-                    elems.append(-1);
-                }
-            }
-
-            elementList_.transfer(elems);
-        }
-
-        // 2. Update faces
-        {
-            DynamicList<label> elems(faceList_.size());
-
-            const labelList& reverseMap = mpm.reverseFaceMap();
-            for (const label facei : faceList_)
-            {
-                if (facei != -1)
-                {
-                    label newFacei = reverseMap[facei];
-                    if (newFacei == -1)
-                    {
-                        // face removed
-                    }
-                    else if (newFacei < -1)
-                    {
-                        // face merged
-                        elems.append(-newFacei - 2);
-                    }
-                    else
-                    {
-                        // valid new face
-                        elems.append(newFacei);
-                    }
-                }
-                else
-                {
-                    // Keep -1 elements
-                    elems.append(-1);
-                }
-            }
-
-            faceList_.transfer(elems);
-        }
-    }
-}
-
-
-void Foam::probes::movePoints(const polyMesh& mesh)
-{
-    DebugInfo<< "probes: movePoints" << endl;
-
-    if (fixedLocations_ && &mesh == &mesh_)
-    {
-        findElements(mesh_);
-    }
-}
-
-
 // ************************************************************************* //

src/sampling/probes/probes.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2022 OpenCFD Ltd.
+    Copyright (C) 2016-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,237 +31,88 @@ Group
     grpUtilitiesFunctionObjects
 
 Description
-    Set of locations to sample.
+    Function object to sample fields at specified internal-mesh locations and
+    write to file.
 
-    Call write() to sample and write files.
-
-    Example of function object specification:
+Usage
+    Minimal example by using \c system/controlDict.functions:
     \verbatim
     probes
     {
+        // Mandatory entries
         type            probes;
         libs            (sampling);
 
-        // Name of the directory for probe data
-        name            probes;
+        fields          (<wordRes>);
+        probeLocations  (<vectorList>);
 
-        // Write at same frequency as fields
-        writeControl    writeTime;
-        writeInterval   1;
+        // Optional entries
+        verbose         <bool>;
+        sampleOnExecute <bool>;
+        fixedLocations  <bool>;
+        includeOutOfBounds <bool>;
+        interpolationScheme <word>;
 
-        // Fields to be probed
-        fields          (U "p.*");
-
-        // Optional: do not recalculate cells if mesh moves
-        fixedLocations  false;
-
-        // Optional: interpolation scheme to use (default is cell)
-        interpolationScheme cellPoint;
-
-        probeLocations
-        (
-            ( 1e-06 0 0.01 )      // at inlet
-            (0.21 -0.20999 0.01)  // at outlet1
-            (0.21 0.20999 0.01)   // at outlet2
-            (0.21 0 0.01)         // at central block
-        );
-
-        // Optional: filter out points that haven't been found. Default
-        //           is to include them (with value -VGREAT)
-        includeOutOfBounds  true;
+        // Inherited entries
+        ...
     }
     \endverbatim
 
-    Entries:
+    where the entries mean:
     \table
-        Property | Description                              | Required | Default
-        type     | Type-name: probes                        | yes |
-        probeLocations | Probe locations                    | yes |
-        fields   | word/regex list of fields to sample      | yes |
-        interpolationScheme | scheme to obtain values       | no  | cell
-        fixedLocations | Do not recalculate cells if mesh moves | no | true
-        includeOutOfBounds | Include out-of-bounds locations | no | true
-        sampleOnExecute | Sample on execution and store results | no | false
+      Property     | Description                           | Type | Reqd | Deflt
+      type         | Type name: probes                     | word |  yes  | -
+      libs         | Library name: sampling                | word |  yes  | -
+      fields       | Names of fields to probe              | wordRes | yes | -
+      probeLocations | Locations of probes              | vectorList | yes | -
+      verbose      | Enable/disable verbose output         | bool | no | false
+      sampleOnExecute | Sample on execution and store results | bool | no <!--
+                   -->                                               | false
+      fixedLocations | Do not recalculate cells if mesh moves | bool | no | true
+      includeOutOfBounds | Include out-of-bounds locations | bool | no | true
+      interpolationScheme | Scheme to obtain values at the points | word <!--
+                   -->                                            | no | cell
     \endtable
 
+    The inherited entries are elaborated in:
+     - \link fvMeshFunctionObject.H \endlink
+     - \link Probes.H \endlink
+     - \link internalProber.H \endlink
+     - \link prober.H \endlink
+
+Note
+  - The \c includeOutOfBounds filters out points that haven't been found.
+    Default is to include them (with value \c -VGREAT).
+
 SourceFiles
     probes.C
-    probesTemplates.C
 
 \*---------------------------------------------------------------------------*/
 
 #ifndef Foam_probes_H
 #define Foam_probes_H
 
-#include "fvMeshFunctionObject.H"
-#include "HashPtrTable.H"
-#include "OFstream.H"
-#include "polyMesh.H"
-#include "pointField.H"
-#include "volFieldsFwd.H"
-#include "surfaceFieldsFwd.H"
-#include "surfaceMesh.H"
-#include "wordRes.H"
-#include "IOobjectList.H"
+#include "Probes.H"
+#include "internalProber.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-// Forward Declarations
-class Time;
-class objectRegistry;
-class dictionary;
-class fvMesh;
-class mapPolyMesh;
-
 /*---------------------------------------------------------------------------*\
                           Class probes Declaration
 \*---------------------------------------------------------------------------*/
 
 class probes
 :
-    public functionObjects::fvMeshFunctionObject,
-    public pointField
+    public Probes<internalProber>
 {
-protected:
+    // Private Data
 
-    // Protected Classes
+        //- Use simpler synonym for the base type
+        using Base = Probes<internalProber>;
 
-        //- Grouping of field names by GeometricField type
-        template<class GeoField>
-        struct fieldGroup : public DynamicList<word> {};
-
-
-    // Data Types
-
-        //- Local control for sampling actions
-        enum sampleActionType : unsigned
-        {
-            ACTION_NONE  = 0,
-            ACTION_WRITE = 0x1,
-            ACTION_STORE = 0x2,
-            ACTION_ALL = 0xF
-        };
-
-
-    // Protected Data
-
-        //- Load fields from files (not from objectRegistry)
-        bool loadFromFiles_;
-
-        //- Fixed locations (default: true)
-        //  Note: set to false for moving mesh calculations where locations
-        //        should move with the mesh
-        bool fixedLocations_;
-
-        //- Include probes that were not found (default: true)
-        bool includeOutOfBounds_;
-
-        //- Output verbosity
-        bool verbose_;
-
-        //- Perform sample actions on execute as well
-        bool onExecute_;
-
-        //- Requested names of fields to probe
-        wordRes fieldSelection_;
-
-        //- Interpolation/sample scheme to obtain values at the points
-        //  Note: only possible when fixedLocations_ is true
-        word samplePointScheme_;
-
-
-    // Calculated
-
-        //- Current list of field names selected for sampling
-        DynamicList<word> selectedFieldNames_;
-
-        //- Categorized scalar/vector/tensor volume fields
-        fieldGroup<volScalarField> scalarFields_;
-        fieldGroup<volVectorField> vectorFields_;
-        fieldGroup<volSphericalTensorField> sphericalTensorFields_;
-        fieldGroup<volSymmTensorField> symmTensorFields_;
-        fieldGroup<volTensorField> tensorFields_;
-
-        //- Categorized scalar/vector/tensor surface fields
-        fieldGroup<surfaceScalarField> surfaceScalarFields_;
-        fieldGroup<surfaceVectorField> surfaceVectorFields_;
-        fieldGroup<surfaceSphericalTensorField> surfaceSphericalTensorFields_;
-        fieldGroup<surfaceSymmTensorField> surfaceSymmTensorFields_;
-        fieldGroup<surfaceTensorField> surfaceTensorFields_;
-
-        //- Cells to be probed (obtained from the locations)
-        labelList elementList_;
-
-        //- Faces to be probed
-        labelList faceList_;
-
-        //- Processor holding the cell or face (-1 if point not found
-        //  on any processor)
-        labelList processor_;
-
-        //- Current open files (non-empty on master only)
-        HashPtrTable<OFstream> probeFilePtrs_;
-
-        //- Patch IDs on which the new probes are located (for patchProbes)
-        labelList patchIDList_;
-
-        //- Original probes location (only used for patchProbes)
-        pointField oldPoints_;
-
-
-    // Protected Member Functions
-
-        //- Find cells and faces containing probes
-        virtual void findElements(const fvMesh& mesh);
-
-        //- Classify field types, close/open file streams
-        //  \return number of fields to sample
-        label prepare(unsigned request);
-
-        //- Get from registry or load from disk
-        template<class GeoField>
-        tmp<GeoField> getOrLoadField(const word& fieldName) const;
-
-        //- Store results: min/max/average/size
-        template<class Type>
-        void storeResults(const word& fieldName, const Field<Type>& values);
-
-
-private:
-
-    // Private Member Functions
-
-        //- Create new streams as required
-        void createProbeFiles(const wordList& fieldNames);
-
-        //- Write field values
-        template<class Type>
-        void writeValues
-        (
-            const word& fieldName,
-            const Field<Type>& values,
-            const scalar timeValue
-        );
-
-        //- Sample and store/write all applicable sampled fields
-        template<class GeoField>
-        void performAction
-        (
-            const fieldGroup<GeoField>& fieldNames,  /* must be sorted */
-            unsigned request
-        );
-
-        //- Perform sampling action with store/write
-        bool performAction(unsigned request);
-
-        //- No copy construct
-        probes(const probes&) = delete;
-
-        //- No copy assignment
-        void operator=(const probes&) = delete;
 
 public:
 
@@ -288,86 +139,19 @@ public:
 
     // Member Functions
 
-        //- Enable/disable verbose output
-        //  \return old value
-        bool verbose(const bool on) noexcept;
+        //- Bring Base::prober into this class's public scope.
+        using Base::prober;
 
-        //- Return names of fields to probe
-        virtual const wordRes& fieldNames() const noexcept
-        {
-            return fieldSelection_;
-        }
-
-        //- Return locations to probe
-        virtual const pointField& probeLocations() const noexcept
-        {
-            return *this;
-        }
-
-        //- Return location for probe i
-        virtual const point& probe(const label i) const
-        {
-            return operator[](i);
-        }
-
-        //- Cells to be probed (obtained from the locations)
-        const labelList& elements() const noexcept
-        {
-            return elementList_;
-        }
-
-        //- Read the probes
+        //- Read the settings from the dictionary
         virtual bool read(const dictionary&);
-
-        //- Sample and store result if the sampleOnExecute is enabled.
-        virtual bool execute();
-
-        //- Sample and write
-        virtual bool write();
-
-        //- Update for changes of mesh
-        virtual void updateMesh(const mapPolyMesh&);
-
-        //- Update for changes of mesh
-        virtual void movePoints(const polyMesh&);
-
-        //- Update for changes of mesh due to readUpdate
-        virtual void readUpdate(const polyMesh::readUpdateState state)
-        {}
-
-
-    // Sampling
-
-        //- Sample a volume field at all locations
-        template<class Type>
-        tmp<Field<Type>> sample(const VolumeField<Type>&) const;
-
-        //- Sample a surface field at all locations
-        template<class Type>
-        tmp<Field<Type>> sample(const SurfaceField<Type>&) const;
-
-        //- Sample a volume field at all locations
-        template<class Type>
-        tmp<Field<Type>> sample(const word& fieldName) const;
-
-        //- Sample a surface field at all locations
-        template<class Type>
-        tmp<Field<Type>> sampleSurfaceField(const word& fieldName) const;
 };
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#ifdef NoRepository
-    #include "probesTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //

src/sampling/probes/probesTemplates.C

@@ -1,274 +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
-    Copyright (C) 2017-2023 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 "probes.H"
-#include "volFields.H"
-#include "surfaceFields.H"
-#include "IOmanip.H"
-#include "interpolation.H"
-#include "SpanStream.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-template<class T>
-struct isNotEqOp
-{
-    void operator()(T& x, const T& y) const
-    {
-        const T unsetVal(-VGREAT*pTraits<T>::one);
-
-        if (x != unsetVal)
-        {
-            // Keep x.
-
-            // Note: should check for y != unsetVal but multiple sample cells
-            // already handled in read().
-        }
-        else
-        {
-            // x is not set. y might be.
-            x = y;
-        }
-    }
-};
-
-} // End namespace Foam
-
-
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-template<class GeoField>
-Foam::tmp<GeoField>
-Foam::probes::getOrLoadField(const word& fieldName) const
-{
-    tmp<GeoField> tfield;
-
-    if (loadFromFiles_)
-    {
-        tfield.emplace
-        (
-            IOobject
-            (
-                fieldName,
-                mesh_.time().timeName(),
-                mesh_.thisDb(),
-                IOobjectOption::MUST_READ,
-                IOobjectOption::NO_WRITE,
-                IOobjectOption::NO_REGISTER
-            ),
-            mesh_
-        );
-    }
-    else
-    {
-        tfield.cref(mesh_.cfindObject<GeoField>(fieldName));
-    }
-
-    return tfield;
-}
-
-
-template<class Type>
-void Foam::probes::storeResults
-(
-    const word& fieldName,
-    const Field<Type>& values
-)
-{
-    const MinMax<Type> limits(values);
-    const Type avgVal = average(values);
-
-    this->setResult("average(" + fieldName + ")", avgVal);
-    this->setResult("min(" + fieldName + ")", limits.min());
-    this->setResult("max(" + fieldName + ")", limits.max());
-    this->setResult("size(" + fieldName + ")", values.size());
-
-    if (verbose_)
-    {
-        Info<< name() << " : " << fieldName << nl
-            << "    avg: " << avgVal << nl
-            << "    min: " << limits.min() << nl
-            << "    max: " << limits.max() << nl << nl;
-    }
-}
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-template<class Type>
-void Foam::probes::writeValues
-(
-    const word& fieldName,
-    const Field<Type>& values,
-    const scalar timeValue
-)
-{
-    if (Pstream::master())
-    {
-        const unsigned int width(IOstream::defaultPrecision() + 7);
-        OFstream& os = *probeFilePtrs_[fieldName];
-
-        os  << setw(width) << timeValue;
-
-        OCharStream buf;
-
-        forAll(values, probei)
-        {
-            if (includeOutOfBounds_ || processor_[probei] != -1)
-            {
-                buf.rewind();
-                buf << values[probei];
-                os  << ' ' << setw(width) << buf.str().data();
-            }
-        }
-        os  << endl;
-    }
-}
-
-
-template<class GeoField>
-void Foam::probes::performAction
-(
-    const fieldGroup<GeoField>& fieldNames,
-    unsigned request
-)
-{
-    for (const word& fieldName : fieldNames)
-    {
-        tmp<GeoField> tfield = getOrLoadField<GeoField>(fieldName);
-
-        if (tfield)
-        {
-            const auto& field = tfield();
-            const scalar timeValue = field.time().timeOutputValue();
-
-            Field<typename GeoField::value_type> values(sample(field));
-
-            this->storeResults(fieldName, values);
-            if (request & ACTION_WRITE)
-            {
-                this->writeValues(fieldName, values, timeValue);
-            }
-        }
-    }
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class Type>
-Foam::tmp<Foam::Field<Type>>
-Foam::probes::sample(const VolumeField<Type>& vField) const
-{
-    const Type unsetVal(-VGREAT*pTraits<Type>::one);
-
-    auto tvalues = tmp<Field<Type>>::New(Field<Type>(this->size(), unsetVal));
-    auto& values = tvalues.ref();
-
-    if (fixedLocations_)
-    {
-        autoPtr<interpolation<Type>> interpPtr
-        (
-            interpolation<Type>::New(samplePointScheme_, vField)
-        );
-
-        forAll(*this, probei)
-        {
-            if (elementList_[probei] >= 0)
-            {
-                const vector& position = operator[](probei);
-
-                values[probei] = interpPtr().interpolate
-                (
-                    position,
-                    elementList_[probei],
-                    -1
-                );
-            }
-        }
-    }
-    else
-    {
-        forAll(*this, probei)
-        {
-            if (elementList_[probei] >= 0)
-            {
-                values[probei] = vField[elementList_[probei]];
-            }
-        }
-    }
-
-    Pstream::listCombineReduce(values, isNotEqOp<Type>());
-
-    return tvalues;
-}
-
-
-template<class Type>
-Foam::tmp<Foam::Field<Type>>
-Foam::probes::sample(const SurfaceField<Type>& sField) const
-{
-    const Type unsetVal(-VGREAT*pTraits<Type>::one);
-
-    auto tvalues = tmp<Field<Type>>::New(Field<Type>(this->size(), unsetVal));
-    auto& values = tvalues.ref();
-
-    forAll(*this, probei)
-    {
-        if (faceList_[probei] >= 0)
-        {
-            values[probei] = sField[faceList_[probei]];
-        }
-    }
-
-    Pstream::listCombineReduce(values, isNotEqOp<Type>());
-
-    return tvalues;
-}
-
-
-template<class Type>
-Foam::tmp<Foam::Field<Type>>
-Foam::probes::sample(const word& fieldName) const
-{
-    return sample(mesh_.lookupObject<VolumeField<Type>>(fieldName));
-}
-
-
-template<class Type>
-Foam::tmp<Foam::Field<Type>>
-Foam::probes::sampleSurfaceField(const word& fieldName) const
-{
-    return sample(mesh_.lookupObject<SurfaceField<Type>>(fieldName));
-}
-
-
-// ************************************************************************* //

src/thermophysicalModels/radiation/radiationModels/radiationModel/radiationModel.H

@@ -243,6 +243,9 @@ public:
             //- Return access to the temperature field
             const volScalarField& T() const noexcept { return T_; }
 
+            //- Return the radiation solver frequency
+            label solverFreq() const noexcept { return solverFreq_; }
+
             //- Source term component (for power of T^4)
             virtual tmp<volScalarField> Rp() const = 0;