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.

applications/test/expressionTemplates/Make/files

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

applications/test/expressionTemplates/Make/options

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

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

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

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

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

src/Allwmake

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

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

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

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

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

src/OpenFOAM/dimensionedTypes/dimensionedType/dimensionedType.H

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

src/OpenFOAM/expressionTemplates/GenericExpression.H

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

src/OpenFOAM/expressionTemplates/boolExpression.H

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

src/OpenFOAM/expressionTemplates/dimensionSetExpression.H

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

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

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

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

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

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

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

src/OpenFOAM/memory/tmp/tmp.H

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

src/OpenFOAM/memory/tmp/tmpI.H

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

src/Pstream/mpi/UPstreamRequest.C

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

src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.C

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C

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

src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H

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

src/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/fused/finiteVolume/fusedGaussGrad.C

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

src/fused/finiteVolume/fusedGaussGrad.H

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

src/fused/finiteVolume/fusedGaussLaplacianScheme.C

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

src/fused/finiteVolume/fusedGaussLaplacianScheme.H

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

src/fused/finiteVolume/fusedGaussLaplacianSchemes.C

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

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

@@ -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/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;
 

tutorials/incompressible/simpleFoam/pitzDaily_fused/Allrun

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

tutorials/incompressible/simpleFoam/pitzDaily_fused/system/controlDict

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

tutorials/incompressible/simpleFoam/pitzDaily_fused/system/fvSolution

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

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

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