TUT: add a multi-region finite-area example tutorial

- the old naming grew organically and was previously restricted
  by avoiding name collisions with volume fields, which resulted
  in some minor consistencies

Common solid shell naming:

  | Keyword  | Default     | Description  | Key(old) | Old Naming     |
  | -------- | ----------- | ------------ | -------- | -------------- |
  | h        | hs (suffix) | thickness    | h        | "h_" + region  |
  | T        | Ts (suffix) | temperature  | -        | "Ts_" + region |

Vibration shells:

  | Keyword  | Default     | Description  | Key(old) | Old Naming     |
  | -------- | ----------- | ------------ | -------- | -------------- |
  |  as      | as (suffix) | acceleration | -        | "as_" + region |
  |  ps      | ps (suffix) | pres on shell| -        | "ps_" + region |
  |  ws      | ws (suffix) | displacement | -        | "ws_" + region |
  |  qs      | qs (suffix  | source field | -        | "qs_" + region |

Other:

  | Keyword  | Default     | Description  | Key(old) |
  | -------- | ----------- | ------------ | -------- |
  | Tprimary | T   | temperature (volume) | T        |

Joule heating: (possibly breaking naming changes)

  | Field    | Name                     | Name(old)                |
  | -------- | ------------------------ | ------------------------ |
  |  V       | \<scope\>:V (suffix)     | \<scope\>:V + region     |
  |  sigma   | \<scope\>:sigma (suffix) | \<scope\>:sigma + region |

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;
-}
-
-
-// ************************************************************************* //

applications/test/faMeshesRegistry/Test-faMeshesRegistry.cxx

@@ -49,8 +49,7 @@ int main(int argc, char *argv[])
 
     Info<< "mesh 0: " << mesh.sortedNames() << nl;
 
-    faMeshesRegistry& reg =
-        const_cast<faMeshesRegistry&>(faMeshesRegistry::New(mesh));
+    auto& reg = const_cast<faMeshesRegistry&>(faMeshesRegistry::New(mesh));
 
     // faMeshesRegistry faReg = faMeshesRegistry(mesh);
 

applications/test/foamToEnsight-check/Make/files

@@ -1,3 +1,3 @@
-foamToEnsight-check.C
+foamToEnsight-check.cxx
 
 EXE = $(FOAM_USER_APPBIN)/foamToEnsight-check

applications/test/foamToEnsight-check/foamToEnsight-check.cxx

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2022-2023 OpenCFD Ltd.
+    Copyright (C) 2022-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -203,6 +203,7 @@ int main(int argc, char *argv[])
     argList::addVerboseOption();
 
     #include "addAllRegionOptions.H"
+    #include "addAllFaRegionOptions.H"
 
     argList::addBoolOption
     (
@@ -251,6 +252,14 @@ int main(int argc, char *argv[])
     // Handle -allRegions, -regions, -region
     #include "getAllRegionOptions.H"
 
+    // Handle -all-area-regions, -area-regions, -area-region
+    #include "getAllFaRegionOptions.H"
+
+    if (!doFiniteArea)
+    {
+        areaRegionNames.clear();  // For consistency
+    }
+
     // ------------------------------------------------------------------------
 
     #include "createNamedMeshes.H"
@@ -259,8 +268,8 @@ int main(int argc, char *argv[])
     /// #include "createMeshAccounting.H"
 
     PtrList<ensightMesh> ensightMeshes(regionNames.size());
-    PtrList<faMesh> meshesFa(regionNames.size());
-    PtrList<ensightFaMesh> ensightMeshesFa(regionNames.size());
+    List<PtrList<faMesh>> meshesFa(regionNames.size());
+    List<PtrList<ensightFaMesh>> ensightMeshesFa(regionNames.size());
 
     forAll(regionNames, regioni)
     {
@@ -273,21 +282,32 @@ int main(int argc, char *argv[])
         );
         ensightMeshes[regioni].verbose(optVerbose);
 
-
-        if (doFiniteArea)
+        if (!doFiniteArea)
         {
-            autoPtr<faMesh> faMeshPtr(faMesh::TryNew(mesh));
+            continue;
+        }
+
+        meshesFa[regioni].resize_null(areaRegionNames.size());
+        ensightMeshesFa[regioni].resize_null(areaRegionNames.size());
+
+        autoPtr<faMesh> faMeshPtr(faMesh::TryNew(mesh));
+
+        forAll(areaRegionNames, areai)
+        {
+            const word& areaName = areaRegionNames[areai];
+
+            autoPtr<faMesh> faMeshPtr(faMesh::TryNew(areaName, mesh));
 
             if (faMeshPtr)
             {
-                meshesFa.set(regioni, std::move(faMeshPtr));
+                meshesFa[regioni].set(areai, std::move(faMeshPtr));
 
-                ensightMeshesFa.set
+                ensightMeshesFa[regioni].set
                 (
-                    regioni,
-                    new ensightFaMesh(meshesFa[regioni])
+                    areai,
+                    new ensightFaMesh(meshesFa[regioni][areai])
                 );
-                ensightMeshesFa[regioni].verbose(optVerbose);
+                ensightMeshesFa[regioni][areai].verbose(optVerbose);
             }
         }
     }
@@ -295,7 +315,7 @@ int main(int argc, char *argv[])
     // ------------------------------------------------------------------------
 
 
-    if (Pstream::master())
+    if (UPstream::master())
     {
         Info<< "Checking " << timeDirs.size() << " time steps" << nl;
     }
@@ -317,17 +337,20 @@ int main(int argc, char *argv[])
             auto& ensMesh = ensightMeshes[regioni];
 
             // Finite-area (can be missing)
-            auto* ensFaMeshPtr = ensightMeshesFa.get(regioni);
+            auto& ensFaMeshes = ensightMeshesFa[regioni];
 
             if (moving)
             {
                 ensMesh.expire();
                 ensMesh.correct();
 
-                if (ensFaMeshPtr)
+                forAll(areaRegionNames, areai)
                 {
-                    ensFaMeshPtr->expire();
-                    ensFaMeshPtr->correct();
+                    if (auto* ensFaMeshPtr = ensFaMeshes.get(areai))
+                    {
+                        ensFaMeshPtr->expire();
+                        ensFaMeshPtr->correct();
+                    }
                 }
             }
 
@@ -340,9 +363,15 @@ int main(int argc, char *argv[])
 
                 printInfo(ensMesh, optVerbose);
 
-                if (ensFaMeshPtr)
+                // finite-area
+                forAll(areaRegionNames, areai)
                 {
-                    printInfo(*ensFaMeshPtr, optVerbose);
+                    auto* ensFaMeshPtr = ensFaMeshes.get(areai);
+
+                    if (ensFaMeshPtr)
+                    {
+                        printInfo(*ensFaMeshPtr, optVerbose);
+                    }
                 }
             }
         }

applications/utilities/finiteArea/checkFaMesh/checkFaMesh.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2021-2023 OpenCFD Ltd.
+    Copyright (C) 2021-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,7 +28,7 @@ Application
     makeFaMesh
 
 Description
-    Check a finiteArea mesh
+    Check a finite-area mesh
 
 Original Authors
     Zeljko Tukovic, FAMENA
@@ -39,6 +39,7 @@ Original Authors
 #include "Time.H"
 #include "argList.H"
 #include "faMesh.H"
+#include "faMeshTools.H"
 #include "polyMesh.H"
 #include "areaFaMesh.H"
 #include "edgeFaMesh.H"
@@ -47,7 +48,7 @@ Original Authors
 #include "processorFaPatch.H"
 #include "foamVtkIndPatchWriter.H"
 #include "foamVtkLineWriter.H"
-#include "faMeshTools.H"
+#include "regionProperties.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -59,7 +60,7 @@ int main(int argc, char *argv[])
 {
     argList::addNote
     (
-        "Check a finiteArea mesh"
+        "Check a finite-area mesh"
     );
 
     argList::addBoolOption
@@ -77,9 +78,15 @@ int main(int argc, char *argv[])
     );
 
     #include "addRegionOption.H"
-    #include "addFaRegionOption.H"
+    #include "addAllFaRegionOptions.H"
     #include "setRootCase.H"
     #include "createTime.H"
+
+    // Handle -all-area-regions, -area-regions, -area-region
+    #include "getAllFaRegionOptions.H"
+
+    // ------------------------------------------------------------------------
+
     #include "createNamedPolyMesh.H"
 
     int geometryOrder(1);
@@ -91,16 +98,41 @@ int main(int argc, char *argv[])
         faMesh::geometryOrder(geometryOrder);
     }
 
-    #include "createNamedFaMesh.H"
-
-    Info<< "Time = " << runTime.timeName() << nl << endl;
-
-    // Mesh information (verbose)
-    faMeshTools::printMeshChecks(aMesh);
-
-    if (args.found("write-vtk"))
+    for (const word& areaName : areaRegionNames)
     {
-        #include "faMeshWriteVTK.H"
+        Info<< "Create faMesh";
+        if (!polyMesh::regionName(areaName).empty())
+        {
+            Info<< " [" << areaName << "]";
+        }
+        Info<< " for time = " << runTime.timeName() << nl;
+
+        autoPtr<faMesh> faMeshPtr(faMesh::TryNew(areaName, mesh));
+
+        if (!faMeshPtr)
+        {
+            Info<< "    ...failed to create area-mesh";
+            if (!polyMesh::regionName(areaName).empty())
+            {
+                Info<< " [" << areaName << "]";
+            }
+            Info<< endl;
+            continue;
+        }
+        else
+        {
+            Info<< endl;
+        }
+
+        const auto& aMesh = faMeshPtr();
+
+        // Mesh information (verbose)
+        faMeshTools::printMeshChecks(aMesh);
+
+        if (args.found("write-vtk"))
+        {
+            #include "faMeshWriteVTK.H"
+        }
     }
 
     Info<< "\nEnd\n" << endl;

applications/utilities/finiteArea/checkFaMesh/faMeshWriteVTK.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021-2023 OpenCFD Ltd.
+    Copyright (C) 2021-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -15,6 +15,16 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
+// The base name for output files
+word vtkBaseFileName(aMesh.regionName());
+if (polyMesh::regionName(vtkBaseFileName).empty())
+{
+    vtkBaseFileName = "finiteArea";
+}
+
+Info<< nl;
+Info<< "Write faMesh in vtk format:" << nl;
+
 {
     // finiteArea - faces
     vtk::uindirectPatchWriter writer
@@ -23,7 +33,7 @@ Description
         // vtk::formatType::INLINE_ASCII,
         fileName
         (
-            aMesh.time().globalPath() / "finiteArea"
+            aMesh.time().globalPath() / vtkBaseFileName
         )
     );
 
@@ -32,7 +42,7 @@ Description
     globalIndex procAddr(aMesh.nFaces());
     labelList cellIDs;
 
-    if (Pstream::master())
+    if (UPstream::master())
     {
         cellIDs.resize(procAddr.totalSize());
         for (const labelRange& range : procAddr.ranges())
@@ -53,8 +63,7 @@ Description
     writer.beginPointData(1);
     writer.write("normal", aMesh.pointAreaNormals());
 
-    Info<< nl
-        << "Wrote faMesh in vtk format: " << writer.output().name() << nl;
+    Info<< "    " << writer.output().name() << nl;
 }
 
 {
@@ -66,7 +75,7 @@ Description
         // vtk::formatType::INLINE_ASCII,
         fileName
         (
-            aMesh.time().globalPath() / "finiteArea-edges"
+            aMesh.time().globalPath() / (vtkBaseFileName + "-edges")
         )
     );
 
@@ -88,8 +97,7 @@ Description
     writer.beginPointData(1);
     writer.write("normal", aMesh.pointAreaNormals());
 
-    Info<< nl
-        << "Wrote faMesh in vtk format: " << writer.output().name() << nl;
+    Info<< "    " << writer.output().name() << nl;
 }
 
 {
@@ -108,7 +116,7 @@ Description
         // vtk::formatType::INLINE_ASCII,
         fileName
         (
-            aMesh.time().globalPath() / "finiteArea-edgesCentres"
+            aMesh.time().globalPath() / (vtkBaseFileName + "-edgesCentres")
         )
     );
 
@@ -134,8 +142,7 @@ Description
         writer.write("normal", fld);
     }
 
-    Info<< nl
-        << "Wrote faMesh in vtk format: " << writer.output().name() << nl;
+    Info<< "    " << writer.output().name() << nl;
 }
 
 

applications/utilities/finiteArea/makeFaMesh/checkPatchTopology.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2023 OpenCFD Ltd.
+    Copyright (C) 2023-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -15,7 +15,7 @@ Description
 
 Input
     mesh (polyMesh)
-    meshDefDict (system/faMeshDefintion)
+    meshDefDict (system/finite-area/faMeshDefinition)
 
 \*---------------------------------------------------------------------------*/
 
@@ -25,13 +25,16 @@ Input
 
     const polyBoundaryMesh& pbm = mesh.boundaryMesh();
 
-    wordRes polyPatchNames;
-    meshDefDict.readIfPresent("polyMeshPatches", polyPatchNames);
+    labelList patchIDs;
 
-    const labelList patchIDs
+    if
     (
-        pbm.indices(polyPatchNames, true)  // useGroups
-    );
+        wordRes select;
+        meshDefDict.readIfPresent("polyMeshPatches", select)
+    )
+    {
+        patchIDs = pbm.indices(select, true);  // useGroups
+    }
 
     label nFaceLabels = 0;
     for (const label patchi : patchIDs)

applications/utilities/finiteArea/makeFaMesh/findMeshDefinitionDict.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021-2024 OpenCFD Ltd.
+    Copyright (C) 2021-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -13,6 +13,18 @@ License
 Description
     Search for the appropriate faMeshDefinition dictionary...
 
+    Expects to find a faMeshDefinition file in a location specified
+    by the command-line option, or one of the standard locations.
+
+    For default area region (region0):
+      (1) system/finite-area/faMeshDefinition
+      (2) system/faMeshDefinition  [legacy location - v2312 and earlier]
+
+    For a named area region:
+      (1) system/finite-area/<area-name>/faMeshDefinition
+      (2) system/finite-area/faMeshDefinition.<area-name>
+          [symmetric with faOptions]
+
 Required Classes
     - Foam::polyMesh
     - Foam::IOdictionary
@@ -24,81 +36,185 @@ Required Variables
     - runTime [Time]
 
 Provided Variables
-    - meshDefDict [IOdictionary]
-    - meshDictPtr [autoPtr<IOdictionary>]
+    - faMeshDefinitionPtr [autoPtr<IOdictionary>]
+
+If the dictionary cannot be found, exits with an error message or
+reports a warning (dry-run mode)
 
 \*---------------------------------------------------------------------------*/
 
 const word dictName("faMeshDefinition");
 
-autoPtr<IOdictionary> meshDictPtr;
+autoPtr<IOdictionary> faMeshDefinitionPtr;
 
 {
-    fileName dictPath;
-    const word& regionDir = Foam::polyMesh::regionName(regionName);
-    const word& areaRegionDir = Foam::polyMesh::regionName(areaRegionName);
+    // Exit unless dry-run?
+    const bool exitOnFailure = (!args.dryRun());
 
-    if (args.readIfPresent("dict", dictPath))
+    const word& regionDir = polyMesh::regionName(regionName);
+    const word& areaRegionDir = polyMesh::regionName(areaRegionName);
+
+    // Canonical location
+    fileName dictPath
+    (
+        runTime.system()/regionDir
+      / faMesh::prefix()/areaRegionDir/dictName
+    );
+
+
+    // Dictionary specified on the command-line ...
+    const bool specified = args.readIfPresent("dict", dictPath);
+
+    if (specified)
     {
-        // Dictionary specified on the command-line ...
-
         if (isDir(dictPath))
         {
             dictPath /= dictName;
         }
     }
-    else if
-    (
-        // Dictionary under system/faMeshDefinition ?
-        // (v2312 and earlier)
 
-        areaRegionDir.empty()
-     && exists
-        (
-            runTime.path()/runTime.caseSystem()
-          / regionDir/faMesh::meshSubDir/dictName
-        )
-    )
-    {
-        // Dictionary present directly in system/ (v2312 and earlier)
-
-        dictPath = runTime.system()/regionDir/dictName;
-    }
-    else
-    {
-        // Use system/finite-area/ directory, with region qualifications
-
-        dictPath =
-        (
-            runTime.system()/regionDir
-          / faMesh::prefix()/areaRegionDir/dictName
-        );
-    }
 
     IOobject meshDictIO
     (
         dictPath,
         runTime,
-        IOobject::MUST_READ,
-        IOobject::NO_WRITE,
-        IOobject::NO_REGISTER,
+        IOobjectOption::MUST_READ,
+        IOobjectOption::NO_WRITE,
+        IOobjectOption::NO_REGISTER,
         true    // is globalObject
     );
 
-    if (!meshDictIO.typeHeaderOk<IOdictionary>(true))
+    refPtr<IOobject> meshDictIOPtr;
+
+    bool foundIOobject = meshDictIO.typeHeaderOk<IOdictionary>(true);
+
+    // For reporting any alternative locations
+    dictPath.clear();
+
+    if (!foundIOobject && !specified)
     {
-        FatalErrorInFunction
-            << meshDictIO.objectPath() << nl
-            << exit(FatalError);
+        if (!areaRegionDir.empty())
+        {
+            // Alternative location
+            // - system/finite-area/faMeshDefinition.<area-name>
+
+            dictPath =
+            (
+                runTime.system()/regionDir
+              / faMesh::prefix()/IOobject::groupName(dictName, areaRegionDir)
+            );
+        }
+        else
+        {
+            // legacy location (v2312 and earlier)
+            // - system/faMeshDefinition
+
+            dictPath = runTime.system()/regionDir/dictName;
+        }
+
+        if (!dictPath.empty())
+        {
+            auto ioptr = autoPtr<IOobject>::New
+            (
+                dictPath,
+                runTime,
+                IOobjectOption::MUST_READ,
+                IOobjectOption::NO_WRITE,
+                IOobjectOption::NO_REGISTER,
+                true    // is globalObject
+            );
+
+            foundIOobject =
+            (
+                ioptr && ioptr->typeHeaderOk<IOdictionary>(true)
+            );
+
+            if (foundIOobject)
+            {
+                // Use if found
+                meshDictIOPtr = std::move(ioptr);
+            }
+
+            // Generally retain dictPath for error messages,
+            // but don't mention the really old legacy location
+            if (areaRegionDir.empty())
+            {
+                dictPath.clear();
+            }
+        }
     }
 
-    Info<< "Creating faMesh from definition: "
-        << meshDictIO.objectRelPath() << endl;
+    // Alternative location or regular location?
+    if (!meshDictIOPtr)
+    {
+        meshDictIOPtr.cref(meshDictIO);
+    }
+    const auto& io = meshDictIOPtr();
 
-    meshDictPtr = autoPtr<IOdictionary>::New(meshDictIO);
+
+    // Handle final success or failure
+    if (foundIOobject)
+    {
+        Info<< "----------------" << nl
+            << "Using area-mesh definition";
+        if (!areaRegionDir.empty())
+        {
+            Info<< " [" << areaRegionName << "]";
+        }
+
+        Info<< nl
+            << "    " << io.objectRelPath() << nl << endl;
+
+        faMeshDefinitionPtr = autoPtr<IOdictionary>::New(io);
+    }
+    else
+    {
+        // Failure
+
+        if (exitOnFailure)
+        {
+            auto& err =
+                FatalErrorInFunction
+                << "Did not find area-mesh definition";
+
+            if (!areaRegionDir.empty())
+            {
+                err<< " [" << areaRegionName << "]";
+            }
+
+            err << nl
+                << "    " << io.objectRelPath() << nl;
+
+            if (!dictPath.empty())
+            {
+                err << "    " << dictPath << nl;
+            }
+
+            FatalError<< exit(FatalError);
+        }
+        else
+        {
+            // dry-run: warning
+            auto& err =
+                Warning << nl
+                << "----------------" << nl
+                << "Did not find area-mesh definition";
+
+            if (!areaRegionDir.empty())
+            {
+                err << " [" << areaRegionName << "]";
+            }
+
+            err << nl
+                << "    " << io.objectRelPath() << nl;
+
+            if (!dictPath.empty())
+            {
+                err << "    " << dictPath << nl;
+            }
+        }
+    }
 }
 
-IOdictionary& meshDefDict = *meshDictPtr;
-
 
 // ************************************************************************* //

applications/utilities/finiteArea/makeFaMesh/makeFaMesh.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2021-2023 OpenCFD Ltd.
+    Copyright (C) 2021-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -40,7 +40,6 @@ Original Authors
 
 #include "Time.H"
 #include "argList.H"
-#include "OSspecific.H"
 #include "faMesh.H"
 #include "faMeshTools.H"
 #include "IOdictionary.H"
@@ -53,6 +52,7 @@ Original Authors
 #include "PtrListOps.H"
 #include "foamVtkLineWriter.H"
 #include "foamVtkIndPatchWriter.H"
+#include "regionProperties.H"
 #include "syncTools.H"
 #include "OBJstream.H"
 
@@ -104,12 +104,14 @@ int main(int argc, char *argv[])
     );
 
     #include "addRegionOption.H"
-    #include "addFaRegionOption.H"
+    #include "addAllFaRegionOptions.H"
+
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createNamedPolyMesh.H"
 
-    #include "getFaRegionOption.H"
+    // Handle -all-area-regions, -area-regions, -area-region
+    #include "getAllFaRegionOptions.H"
 
     const bool doDecompose = !args.found("no-decompose");
     const bool doDecompFields = !args.found("no-fields");
@@ -123,56 +125,83 @@ int main(int argc, char *argv[])
         Info<< "Skip decompose of finiteArea fields" << nl;
     }
 
-    // Reading faMeshDefinition dictionary
-    #include "findMeshDefinitionDict.H"
+    // ------------------------------------------------------------------------
 
-    // Inject/overwrite name for optional 'empty' patch
-    word patchName;
-    if (args.readIfPresent("empty-patch", patchName))
+    for (const word& areaRegionName : areaRegionNames)
     {
-        meshDefDict.add("emptyPatch", patchName, true);
-    }
+        // Reading faMeshDefinition dictionary
+        #include "findMeshDefinitionDict.H"
 
-    // Preliminary checks
-    #include "checkPatchTopology.H"
+        if (!faMeshDefinitionPtr)
+        {
+            if (args.dryRun())
+            {
+                break;
+            }
+            else
+            {
+                FatalErrorInFunction
+                    << "Did not find area-mesh definition"<< nl
+                    << exit(FatalError);
+            }
+        }
 
-    Info << "Create areaMesh";
-    if (!Foam::polyMesh::regionName(areaRegionName).empty())
-    {
-        Foam::Info << ' ' << areaRegionName;
-    }
-    Info << " for polyMesh at time = " << runTime.timeName() << nl;
+        auto& meshDefDict = (*faMeshDefinitionPtr);
 
-    // Create
-    faMesh aMesh(areaRegionName, mesh, meshDefDict);
 
-    // Mesh information (less verbose)
-    faMeshTools::printMeshChecks(aMesh, 0);
+        // Inject/overwrite name for optional 'empty' patch
+        if (word name; args.readIfPresent("empty-patch", name))
+        {
+            meshDefDict.add("emptyPatch", name, true);
+        }
 
-    if (args.found("write-edges-obj"))
-    {
-        #include "faMeshWriteEdgesOBJ.H"
-    }
+        // Preliminary checks
+        #include "checkPatchTopology.H"
 
-    if (args.found("write-vtk"))
-    {
-        #include "faMeshWriteVTK.H"
-    }
+        Info<< "Create area-mesh";
+        if (!polyMesh::regionName(areaRegionName).empty())
+        {
+            Info<< " [" << areaRegionName << "]";
+        }
+        Info<< " with polyMesh at time = " << runTime.timeName() << nl;
 
-    if (args.dryRun())
-    {
-        Info<< "\ndry-run: not writing mesh or decomposing fields\n" << nl;
-    }
-    else
-    {
-        // More precision (for points data)
-        IOstream::minPrecision(10);
+        // Create
+        faMesh aMesh(areaRegionName, mesh, meshDefDict);
 
-        Info<< nl << "Write finite area mesh." << nl;
-        aMesh.write();
+        // Mesh information (less verbose)
+        faMeshTools::printMeshChecks(aMesh, 0);
 
-        Info<< endl;
-        #include "decomposeFaFields.H"
+        if (args.found("write-edges-obj"))
+        {
+            #include "faMeshWriteEdgesOBJ.H"
+        }
+
+        if (args.found("write-vtk"))
+        {
+            #include "faMeshWriteVTK.H"
+        }
+
+        if (args.dryRun())
+        {
+            Info<< "\ndry-run: not writing mesh or decomposing fields\n" << nl;
+        }
+        else
+        {
+            // More precision (for points data)
+            IOstream::minPrecision(10);
+
+            Info<< nl << "Write finite area mesh";
+            if (const word& name = aMesh.regionName(); !name.empty())
+            {
+                Info<< " [" << name << "]";
+            }
+            Info<< nl;
+
+            aMesh.write();
+
+            Info<< endl;
+            #include "decomposeFaFields.H"
+        }
     }
 
     Info<< "\nEnd\n" << endl;

applications/utilities/parallelProcessing/decomposePar/decomposePar.C

@@ -334,6 +334,7 @@ int main(int argc, char *argv[])
     );
 
     #include "addAllRegionOptions.H"
+    #include "addAllFaRegionOptions.H"
 
     argList::addDryRunOption
     (
@@ -442,7 +443,6 @@ int main(int argc, char *argv[])
     bool decomposeFieldsOnly = args.found("fields");
     bool forceOverwrite      = args.found("force");
 
-
     // Set time from database
     #include "createTime.H"
 
@@ -491,9 +491,17 @@ int main(int argc, char *argv[])
         decompDictFile = runTime.globalPath()/decompDictFile;
     }
 
-    // Get region names
+    // Handle -allRegions, -regions, -region
     #include "getAllRegionOptions.H"
 
+    // Handle -all-area-regions, -area-regions, -area-region
+    #include "getAllFaRegionOptions.H"
+
+    if (!doFiniteArea)
+    {
+        areaRegionNames.clear();  // For consistency
+    }
+
     const bool optRegions =
         (regionNames.size() != 1 || regionNames[0] != polyMesh::defaultRegion);
 
@@ -813,56 +821,76 @@ int main(int argc, char *argv[])
 
                 // Field objects at this time
                 IOobjectList objects;
-                IOobjectList faObjects;
+
+                // faMesh fields - can have multiple finite-area per volume
+                HashTable<IOobjectList> faObjects;
 
                 if (doDecompFields)
                 {
                     // List of volume mesh objects for this instance
                     objects = IOobjectList(mesh, runTime.timeName());
 
-                    // List of area mesh objects (assuming single region)
-                    faObjects = IOobjectList
-                    (
-                        mesh.time(),
-                        runTime.timeName(),
-                        faMesh::dbDir(mesh, word::null),
-                        IOobjectOption::NO_REGISTER
-                    );
-
                     // Ignore generated fields: (cellDist)
                     objects.remove("cellDist");
+
+                    // Lists of finite-area fields
+                    faObjects.reserve(areaRegionNames.size());
+
+                    for (const word& areaName : areaRegionNames)
+                    {
+                        // The finite-area objects for this area region
+                        IOobjectList objs
+                        (
+                            faMesh::Registry(mesh),
+                            runTime.timeName(),
+                            polyMesh::regionName(areaName),
+                            IOobjectOption::NO_REGISTER
+                        );
+
+                        if (!objs.empty())
+                        {
+                            faObjects.emplace_set(areaName, std::move(objs));
+                        }
+                    }
                 }
 
                 // Finite area handling
-                autoPtr<faMeshDecomposition> faMeshDecompPtr;
+                // - all area regions use the same volume decomposition
+
+                HashPtrTable<faMeshDecomposition> faMeshDecompHashes;
                 if (doFiniteArea)
                 {
                     const word boundaryInst =
                         mesh.time().findInstance(mesh.meshDir(), "boundary");
 
-                    IOobject io
-                    (
-                        "faBoundary",
-                        boundaryInst,
-                        faMesh::meshDir(mesh, word::null),
-                        mesh.time(),
-                        IOobject::READ_IF_PRESENT,
-                        IOobject::NO_WRITE,
-                        IOobject::NO_REGISTER
-                    );
-
-                    if (io.typeHeaderOk<faBoundaryMesh>(true))
+                    for (const word& areaName : areaRegionNames)
                     {
-                        // Always based on the volume decomposition!
-                        faMeshDecompPtr.reset
+                        IOobject io
                         (
-                            new faMeshDecomposition
-                            (
-                                mesh,
-                                mesh.nProcs(),
-                                mesh.model()
-                            )
+                            "faBoundary",
+                            boundaryInst,
+                            faMesh::meshDir(mesh, areaName),
+                            mesh.time(),
+                            IOobject::READ_IF_PRESENT,
+                            IOobject::NO_WRITE,
+                            IOobject::NO_REGISTER
                         );
+
+                        if (io.typeHeaderOk<faBoundaryMesh>(true))
+                        {
+                            // Always based on the volume decomposition!
+                            faMeshDecompHashes.set
+                            (
+                                areaName,
+                                autoPtr<faMeshDecomposition>::New
+                                (
+                                    areaName,
+                                    mesh,
+                                    mesh.nProcs(),
+                                    mesh.model()
+                                )
+                            );
+                        }
                     }
                 }
 
@@ -1066,7 +1094,7 @@ int main(int argc, char *argv[])
 
                     processorDb.setTime(runTime);
 
-                    // read the mesh
+                    // Read the mesh
                     if (!procMeshList.set(proci))
                     {
                         procMeshList.set
@@ -1273,12 +1301,15 @@ int main(int argc, char *argv[])
                 }
 
 
-                // Finite area mesh and field decomposition
-                if (faMeshDecompPtr)
+                // Finite-area mesh and field decomposition
+                for (auto& iter : faMeshDecompHashes.sorted())
                 {
-                    Info<< "\nFinite area mesh decomposition" << endl;
+                    const word& areaName = iter.key();
 
-                    faMeshDecomposition& aMesh = faMeshDecompPtr();
+                    faMeshDecomposition& aMesh = *(iter.val());
+
+                    Info<< "\nFinite area mesh decomposition: "
+                        << areaName << endl;
 
                     aMesh.decomposeMesh();
                     aMesh.writeDecomposition();
@@ -1289,9 +1320,13 @@ int main(int argc, char *argv[])
 
                     faFieldDecomposer::fieldsCache areaFieldCache;
 
-                    if (doDecompFields)
+                    if
+                    (
+                        const auto objs = faObjects.cfind(areaName);
+                        doDecompFields && objs.good()
+                    )
                     {
-                        areaFieldCache.readAllFields(aMesh, faObjects);
+                        areaFieldCache.readAllFields(aMesh, objs.val());
                     }
 
                     const label nAreaFields = areaFieldCache.size();
@@ -1322,7 +1357,7 @@ int main(int argc, char *argv[])
 
                         processorDb.setTime(runTime);
 
-                        // Read the mesh
+                        // Read the volume mesh
                         fvMesh procFvMesh
                         (
                             IOobject
@@ -1333,7 +1368,7 @@ int main(int argc, char *argv[])
                             )
                         );
 
-                        faMesh procMesh(procFvMesh);
+                        faMesh procMesh(areaName, procFvMesh);
 
                         // // Does not work.  HJ, 15/Aug/2017
                         // const labelIOList& faceProcAddressing =
@@ -1392,7 +1427,11 @@ int main(int argc, char *argv[])
                             boundaryProcAddressing
                         );
 
-                        areaFieldCache.decomposeAllFields(fieldDecomposer);
+                        areaFieldCache.decomposeAllFields
+                        (
+                            fieldDecomposer,
+                            args.verbose()  // report
+                        );
                     }
                 }
             }

applications/utilities/parallelProcessing/reconstructPar/reconstructPar.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.
@@ -91,6 +91,7 @@ int main(int argc, char *argv[])
     argList::noParallel();
 
     #include "addAllRegionOptions.H"
+    #include "addAllFaRegionOptions.H"
 
     argList::addVerboseOption();
     argList::addOption
@@ -201,9 +202,17 @@ int main(int argc, char *argv[])
 
     const bool newTimes = args.found("newTimes");
 
-    // Get region names
+    // Handle -allRegions, -regions, -region
     #include "getAllRegionOptions.H"
 
+    // Handle -all-area-regions, -area-regions, -area-region
+    #include "getAllFaRegionOptions.H"
+
+    if (!doFiniteArea)
+    {
+        areaRegionNames.clear();  // For consistency
+    }
+
     // Determine the processor count
     label nProcs{0};
 
@@ -391,19 +400,30 @@ int main(int argc, char *argv[])
                 IOobjectOption::NO_REGISTER
             );
 
-            IOobjectList faObjects;
+            // The finite-area fields (multiple finite-area per volume)
+            HashTable<IOobjectList> faObjects;
 
             if (doFiniteArea && doFields)
             {
-                // List of area mesh objects (assuming single region)
-                // - scan on processor0
-                faObjects = IOobjectList
-                (
-                    procMeshes.meshes()[0],
-                    databases[0].timeName(),
-                    faMesh::dbDir(word::null),  // local relative to mesh
-                    IOobjectOption::NO_REGISTER
-                );
+                // Lists of finite-area fields - scan on processor0
+                faObjects.reserve(areaRegionNames.size());
+
+                for (const word& areaName : areaRegionNames)
+                {
+                    // The finite-area objects for this area region
+                    IOobjectList objs
+                    (
+                        procMeshes.meshes()[0],
+                        databases[0].timeName(),
+                        faMesh::dbDir(areaName),  // local relative to mesh
+                        IOobjectOption::NO_REGISTER
+                    );
+
+                    if (!objs.empty())
+                    {
+                        faObjects.emplace_set(areaName, std::move(objs));
+                    }
+                }
             }
 
             if (doFields)
@@ -553,42 +573,60 @@ int main(int argc, char *argv[])
             }
 
 
-            // If there are any FA fields, reconstruct them
-
-            if (!doFiniteArea)
+            // Reconstruct any finite-area fields
+            if (doFiniteArea)
             {
+                bool hadFaFields = false;
+                for (const word& areaName : areaRegionNames)
+                {
+                    const auto objs = faObjects.cfind(areaName);
+                    if (!objs.good())
+                    {
+                        continue;
+                    }
+                    const auto& faObjs = objs.val();
+
+                    if
+                    (
+                        !faObjs.count(fieldTypes::is_area)
+                     && !faObjs.count<edgeScalarField>()
+                    )
+                    {
+                        continue;
+                    }
+
+                    hadFaFields = true;
+
+                    Info<< "Reconstructing finite-area fields ["
+                        << polyMesh::regionName(areaName)
+                        << "]" << nl << endl;
+
+                    const faMesh aMesh(areaName, mesh);
+
+                    processorFaMeshes procFaMeshes
+                    (
+                        procMeshes.meshes(),
+                        areaName
+                    );
+
+                    faFieldReconstructor reconstructor
+                    (
+                        aMesh,
+                        procFaMeshes.meshes(),
+                        procFaMeshes.edgeProcAddressing(),
+                        procFaMeshes.faceProcAddressing(),
+                        procFaMeshes.boundaryProcAddressing()
+                    );
+
+                    reconstructor.reconstructAllFields(faObjs);
+                }
+
+                if (!hadFaFields)
+                {
+                    Info << "No finite-area fields" << nl << endl;
+                }
             }
-            else if
-            (
-                faObjects.count<areaScalarField>()
-             || faObjects.count<areaVectorField>()
-             || faObjects.count<areaSphericalTensorField>()
-             || faObjects.count<areaSymmTensorField>()
-             || faObjects.count<areaTensorField>()
-             || faObjects.count<edgeScalarField>()
-            )
-            {
-                Info << "Reconstructing FA fields" << nl << endl;
 
-                faMesh aMesh(mesh);
-
-                processorFaMeshes procFaMeshes(procMeshes.meshes());
-
-                faFieldReconstructor reconstructor
-                (
-                    aMesh,
-                    procFaMeshes.meshes(),
-                    procFaMeshes.edgeProcAddressing(),
-                    procFaMeshes.faceProcAddressing(),
-                    procFaMeshes.boundaryProcAddressing()
-                );
-
-                reconstructor.reconstructAllFields(faObjects);
-            }
-            else
-            {
-                Info << "No FA fields" << nl << endl;
-            }
 
             if (doReconstructSets)
             {

applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C

@@ -119,17 +119,10 @@ autoPtr<faceCoupleInfo> determineCoupledFaces
         const polyBoundaryMesh& masterPatches = masterMesh.boundaryMesh();
 
 
-        DynamicList<label> masterFaces
-        (
-            masterMesh.nFaces()
-          - masterMesh.nInternalFaces()
-        );
+        DynamicList<label> masterFaces(masterMesh.nBoundaryFaces());
 
-
-        forAll(masterPatches, patchi)
+        for (const polyPatch& pp : masterPatches)
         {
-            const polyPatch& pp = masterPatches[patchi];
-
             if (isA<processorPolyPatch>(pp))
             {
                 for
@@ -139,11 +132,8 @@ autoPtr<faceCoupleInfo> determineCoupledFaces
                     proci++
                 )
                 {
-                    const string toProcString("to" + name(proci));
-                    if (
-                        pp.name().rfind(toProcString)
-                     == (pp.name().size()-toProcString.size())
-                    )
+                    const string toProcString("to" + Foam::name(proci));
+                    if (pp.name().ends_with(toProcString))
                     {
                         label meshFacei = pp.start();
                         forAll(pp, i)
@@ -156,7 +146,6 @@ autoPtr<faceCoupleInfo> determineCoupledFaces
 
             }
         }
-        masterFaces.shrink();
 
 
         // Pick up all patches on meshToAdd ending in "procBoundaryDDDtoYYY"
@@ -164,16 +153,10 @@ autoPtr<faceCoupleInfo> determineCoupledFaces
 
         const polyBoundaryMesh& addPatches = meshToAdd.boundaryMesh();
 
-        DynamicList<label> addFaces
-        (
-            meshToAdd.nFaces()
-          - meshToAdd.nInternalFaces()
-        );
+        DynamicList<label> addFaces(meshToAdd.nBoundaryFaces());
 
-        forAll(addPatches, patchi)
+        for (const polyPatch& pp : addPatches)
         {
-            const polyPatch& pp = addPatches[patchi];
-
             if (isA<processorPolyPatch>(pp))
             {
                 bool isConnected = false;
@@ -215,7 +198,6 @@ autoPtr<faceCoupleInfo> determineCoupledFaces
                 }
             }
         }
-        addFaces.shrink();
 
         return autoPtr<faceCoupleInfo>::New
         (
@@ -502,15 +484,14 @@ void writeMaps
 
     Info<< "    pointProcAddressing" << endl;
     ioAddr.rename("pointProcAddressing");
-    IOList<label>::writeContents(ioAddr, pointProcAddressing);
+    labelIOList::writeContents(ioAddr, pointProcAddressing);
 
     // From processor face to reconstructed mesh face
-    Info<< "    faceProcAddressing" << endl;
-    ioAddr.rename("faceProcAddressing");
-    labelIOList faceProcAddr(ioAddr, faceProcAddressing);
+    // - add turning index to faceProcAddressing.
+    //   See reconstructPar for meaning of turning index.
+
+    labelList faceProcAddr(faceProcAddressing);
 
-    // Now add turning index to faceProcAddressing.
-    // See reconstructPar for meaning of turning index.
     forAll(faceProcAddr, procFacei)
     {
         const label masterFacei = faceProcAddr[procFacei];
@@ -545,19 +526,21 @@ void writeMaps
         }
     }
 
-    faceProcAddr.write();
-
+    Info<< "    faceProcAddressing" << endl;
+    ioAddr.rename("faceProcAddressing");
+    labelIOList::writeContents(ioAddr, faceProcAddr);
+    faceProcAddr.clear();
 
     // From processor cell to reconstructed mesh cell
     Info<< "    cellProcAddressing" << endl;
     ioAddr.rename("cellProcAddressing");
-    IOList<label>::writeContents(ioAddr, cellProcAddressing);
+    labelIOList::writeContents(ioAddr, cellProcAddressing);
 
 
     // From processor patch to reconstructed mesh patch
     Info<< "    boundaryProcAddressing" << endl;
     ioAddr.rename("boundaryProcAddressing");
-    IOList<label>::writeContents(ioAddr, boundProcAddressing);
+    labelIOList::writeContents(ioAddr, boundProcAddressing);
 
     Info<< endl;
 }
@@ -645,7 +628,8 @@ void sortFaEdgeMapping
 {
     // From faMeshReconstructor.C - edge shuffling on patches
 
-    Map<label> remapGlobal(2*onePatch.nEdges());
+    Map<label> remapGlobal;
+    remapGlobal.reserve(onePatch.nEdges());
     for (label edgei = 0; edgei < onePatch.nInternalEdges(); ++edgei)
     {
         remapGlobal.insert(edgei, remapGlobal.size());
@@ -779,6 +763,11 @@ int main(int argc, char *argv[])
     );
 
     #include "addAllRegionOptions.H"
+    #include "addAllFaRegionOptions.H"
+
+    // Prevent volume fields [with regionFaModels] from incidental
+    // triggering finite-area
+    regionModels::allowFaModels(false);
 
     // Prevent volume BCs from triggering finite-area
     regionModels::allowFaModels(false);
@@ -841,9 +830,17 @@ int main(int argc, char *argv[])
         }
     }
 
-    // Get region names
+    // Handle -allRegions, -regions, -region
     #include "getAllRegionOptions.H"
 
+    // Handle -all-area-regions, -area-regions, -area-region
+    #include "getAllFaRegionOptions.H"
+
+    if (!doFiniteArea)
+    {
+        areaRegionNames.clear();  // For consistency
+    }
+
     // Determine the processor count
     label nProcs{0};
 
@@ -935,7 +932,8 @@ int main(int argc, char *argv[])
                 polyMesh::meshDir(regionName),
                 databases[0],
                 IOobject::NO_READ,
-                IOobject::NO_WRITE
+                IOobject::NO_WRITE,
+                IOobject::NO_REGISTER
             );
 
             // Problem: faceCompactIOList recognises both 'faceList' and
@@ -1405,8 +1403,12 @@ int main(int argc, char *argv[])
 
 
             // Read finite-area
-            PtrList<faMesh> procFaMeshes(databases.size());
+            // For each named area region that exists create a HashTable
+            // entry that will contain the PtrList for all processors
+
             PtrList<polyMesh> procMeshes(databases.size());
+            HashTable<PtrList<faMesh>> procAreaRegionMeshes;
+            procAreaRegionMeshes.reserve(areaRegionNames.size());
 
             forAll(databases, proci)
             {
@@ -1414,20 +1416,16 @@ int main(int argc, char *argv[])
                     << "Read processor mesh: "
                     << (databases[proci].caseName()/regionDir) << endl;
 
-                procMeshes.set
+                const polyMesh& procMesh = procMeshes.emplace
                 (
                     proci,
-                    new polyMesh
+                    IOobject
                     (
-                        IOobject
-                        (
-                            regionName,
-                            databases[proci].timeName(),
-                            databases[proci]
-                        )
+                        regionName,
+                        databases[proci].timeName(),
+                        databases[proci]
                     )
                 );
-                const polyMesh& procMesh = procMeshes[proci];
 
                 writeMaps
                 (
@@ -1450,38 +1448,67 @@ int main(int argc, char *argv[])
                             "boundary"
                         );
 
-                    IOobject io
-                    (
-                        "faBoundary",
-                        boundaryInst,
-                        faMesh::meshDir(procMesh, word::null),
-                        procMesh.time(),
-                        IOobject::READ_IF_PRESENT,
-                        IOobject::NO_WRITE,
-                        IOobject::NO_REGISTER
-                    );
-
-                    if (io.typeHeaderOk<faBoundaryMesh>(true))
+                    for (const word& areaName : areaRegionNames)
                     {
-                        // Always based on the volume decomposition!
-                        procFaMeshes.set(proci, new faMesh(procMesh));
+                        IOobject io
+                        (
+                            "faBoundary",
+                            boundaryInst,
+                            faMesh::meshDir(procMesh, areaName),
+                            procMesh.time(),
+                            IOobject::READ_IF_PRESENT,
+                            IOobject::NO_WRITE,
+                            IOobject::NO_REGISTER
+                        );
+
+                        if (io.typeHeaderOk<faBoundaryMesh>(true))
+                        {
+                            // Always based on the volume decomposition!
+                            auto& procFaMeshes = procAreaRegionMeshes(areaName);
+                            procFaMeshes.resize(databases.size());
+
+                            procFaMeshes.set
+                            (
+                                proci,
+                                new faMesh(areaName, procMesh)
+                            );
+                        }
                     }
                 }
             }
 
 
-            // Finite-area mapping
-            doFiniteArea = false;
-            forAll(procFaMeshes, proci)
+            // A finite-area mapping exists if procFaMeshes was filled
+
+            // Re-read reconstructed polyMesh. Note: could probably be avoided
+            // by merging into loops above.
+            refPtr<polyMesh> masterPolyMeshPtr;
+
+            if (!procAreaRegionMeshes.empty())
             {
-                if (procFaMeshes.set(proci))
-                {
-                    doFiniteArea = true;
-                }
+                masterPolyMeshPtr.reset
+                (
+                    new polyMesh
+                    (
+                        IOobject
+                        (
+                            regionName,
+                            runTime.timeName(),
+                            runTime
+                        ),
+                        true
+                    )
+                );
             }
 
-            if (doFiniteArea)
+            // Process any finite-area meshes
+            for (const auto& iter : procAreaRegionMeshes.csorted())
             {
+                const auto& areaName = iter.key();
+                const auto& procFaMeshes = iter.val();
+
+                const polyMesh& masterMesh = masterPolyMeshPtr();
+
                 // Addressing from processor to reconstructed case
                 labelListList faFaceProcAddressing(nProcs);
                 labelListList faEdgeProcAddressing(nProcs);
@@ -1499,32 +1526,10 @@ int main(int argc, char *argv[])
 
                     faBoundProcAddressing[proci] = identity(bm.size());
                     // Mark processor patches
-                    for
-                    (
-                        label patchi = bm.nNonProcessor();
-                        patchi < bm.size();
-                        ++patchi
-                    )
-                    {
-                        faBoundProcAddressing[proci][patchi] = -1;
-                    }
+                    faBoundProcAddressing[proci].slice(bm.nNonProcessor()) = -1;
                 }
 
 
-                // Re-read reconstructed polyMesh. Note: could probably be avoided
-                // by merging into loops above.
-                const polyMesh masterMesh
-                (
-                    IOobject
-                    (
-                        regionName,
-                        runTime.timeName(),
-                        runTime
-                    ),
-                    true
-                );
-
-
                 // faceProcAddressing
                 // ~~~~~~~~~~~~~~~~~~
 
@@ -1559,13 +1564,13 @@ int main(int argc, char *argv[])
                 // Construct without patches
                 faMesh masterFaMesh
                 (
+                    areaName,
                     masterMesh,
                     std::move(masterFaceLabels),
                     io
                 );
 
-                const uindirectPrimitivePatch& masterPatch =
-                    masterFaMesh.patch();
+                const auto& masterPatch = masterFaMesh.patch();
 
 
                 // pointProcAddressing
@@ -1577,7 +1582,7 @@ int main(int argc, char *argv[])
                     const auto& procPatch = procFaMeshes[proci].patch();
                     const auto& mp = procPatch.meshPoints();
 
-                    labelList& pointAddr = faPointProcAddressing[proci];
+                    auto& pointAddr = faPointProcAddressing[proci];
                     pointAddr.resize_nocopy(mp.size());
 
                     forAll(mp, i)
@@ -1626,8 +1631,7 @@ int main(int argc, char *argv[])
                 label nPatches = 0;
                 forAll(completePatches, patchi)
                 {
-                    const labelList& patchEdgeLabels =
-                        singlePatchEdgeLabels[patchi];
+                    const auto& patchEdgeLabels = singlePatchEdgeLabels[patchi];
 
                     const faPatch& fap = procMesh0.boundary()[patchi];
 
@@ -1658,11 +1662,11 @@ int main(int argc, char *argv[])
 
                 // Serial mesh - no parallel communication
 
-                const bool oldParRun = Pstream::parRun(false);
+                const bool oldParRun = UPstream::parRun(false);
 
                 masterFaMesh.addFaPatches(completePatches);
 
-                Pstream::parRun(oldParRun);  // Restore parallel state
+                UPstream::parRun(oldParRun);  // Restore parallel state
 
 
                 // Write mesh & individual addressing

applications/utilities/parallelProcessing/redistributePar/loadOrCreateMesh.C

@@ -33,99 +33,266 @@ License
 #include "decomposedBlockData.H"
 #include "IFstream.H"
 
-// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
 
-bool Foam::checkFileExistence(const fileName& fName)
+// Trimmed-down version of lookupAndCacheProcessorsPath
+// with Foam::exists() check. No caching.
+
+// Check for two conditions:
+// - file has to exist
+// - if collated the entry has to exist inside the file
+
+// Note: bypass fileOperation::filePath(IOobject&) since has problems
+//       with going to a different number of processors
+//       (in collated format). Use file-based searching instead
+
+namespace Foam
 {
-    // Trimmed-down version of lookupAndCacheProcessorsPath
-    // with Foam::exists() check. No caching.
 
-    // Check for two conditions:
-    // - file has to exist
-    // - if collated the entry has to exist inside the file
-
-    // Note: bypass fileOperation::filePath(IOobject&) since has problems
-    //       with going to a different number of processors
-    //       (in collated format). Use file-based searching instead
-
-    const auto& handler = Foam::fileHandler();
-    typedef fileOperation::procRangeType procRangeType;
-
-    fileName path, pDir, local;
-    procRangeType group;
-    label numProcs;
-    const label proci =
-        fileOperation::splitProcessorPath
-        (fName, path, pDir, local, group, numProcs);
+// If indeed collated format:
+// Collect block-number in individual filenames
+// (might differ on different processors)
+static bool checkFileExistenceCollated
+(
+    const Foam::fileOperation& handler,
+    const Foam::fileName& fName
+)
+{
+    // using namespace Foam;
 
     bool found = false;
-
-    if (proci != -1)
     {
-        // Read all directories to see any beginning with processor
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        const label handlerComm = handler.comm();
 
-        const fileNameList dirEntries
-        (
-            handler.readDir(path, fileName::Type::DIRECTORY)
-        );
+        const label globalProci = UPstream::myProcNo(UPstream::worldComm);
+        const label handlerProci = UPstream::myProcNo(handlerComm);
+        const label nHandlerProcs = UPstream::nProcs(handlerComm);
 
-        // Extract info from processorN or processorsNN
-        // - highest processor number
-        // - directory+offset containing data for proci
-
-        // label nProcs = 0;
-        for (const fileName& dirN : dirEntries)
+        // Determine my local block number
+        label myBlockNumber = -1;
         {
-            // Analyse directory name
-            label rNum(-1);
-            const label readProci =
-                fileOperation::detectProcessorPath(dirN, group, &rNum);
+            fileOperation::procRangeType group;
+            label proci = fileOperation::detectProcessorPath(fName, group);
 
-            if (proci == readProci)
+            if (proci == -1 && group.empty())
             {
-                // Found "processorN"
-                if (Foam::exists(path/dirN/local))
-                {
-                    found = true;
-                    break;
-                }
+                // 'processorsXXX' format so contains all ranks
+                // according to worldComm
+                myBlockNumber = globalProci;
             }
-            else if (rNum != -1)
+            else
             {
-                // "processorsNN" or "processorsNN_start-end"
-                if (group.empty())
-                {
-                    // "processorsNN"
-                    if (proci < rNum && Foam::exists(path/dirN/local))
-                    {
-                        found = true;
-                        break;
-                    }
-                }
-                else if (group.contains(proci))
-                {
-                    // "processorsNN_start-end"
-                    // - save the local proc offset
-
-                    if (Foam::exists(path/dirN/local))
-                    {
-                        found = true;
-                        break;
-                    }
-                }
+                // 'processorsXXX_n-m' format so check for relative rank
+                myBlockNumber = handlerProci;
             }
         }
-    }
 
-    if (!found)
-    {
-        found = Foam::exists(fName);
+        // Since we are streaming anyhow, could also pack as tuple:
+        // Tuple2<fileName, label>
+
+
+        // Collect file names on master of local communicator
+        const fileNameList fNames
+        (
+            Pstream::listGatherValues
+            (
+                fName,
+                handlerComm,
+                UPstream::msgType()
+            )
+        );
+
+        // Collect block numbers on master of local communicator
+        const labelList myBlockNumbers
+        (
+            Pstream::listGatherValues
+            (
+                myBlockNumber,
+                handlerComm,
+                UPstream::msgType()
+            )
+        );
+
+        // Determine for all whether the filename exists in the collated file.
+        boolList allFound;
+
+        if (UPstream::master(handlerComm))
+        {
+            allFound.resize(nHandlerProcs, false);
+
+            // Store nBlocks and index of file that was used for nBlocks
+            label nBlocks = -1;
+            label blockRanki = -1;
+            forAll(fNames, ranki)
+            {
+                if
+                (
+                    blockRanki == -1
+                 || (fNames[ranki] != fNames[blockRanki])
+                )
+                {
+                    blockRanki = ranki;
+                    IFstream is(fNames[ranki]);
+                    nBlocks = decomposedBlockData::getNumBlocks(is);
+                }
+
+                allFound[ranki] = (myBlockNumbers[ranki] < nBlocks);
+            }
+        }
+
+        // Scatter using the handler communicator
+        found = Pstream::listScatterValues
+        (
+            allFound,
+            handlerComm,
+            UPstream::msgType()
+        );
     }
 
     return found;
 }
 
+} // End namespace
+
+
+// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
+
+Foam::bitSet Foam::haveProcessorFile
+(
+    const word& name,           // eg "faces"
+    const fileName& instance,   // eg "constant"
+    const fileName& local,      // eg, polyMesh
+    const Time& runTime,
+    const bool verbose
+)
+{
+    const auto& handler = Foam::fileHandler();
+
+    const fileName fName
+    (
+        handler.filePath(runTime.path()/instance/local/name)
+    );
+
+    bool found = handler.isFile(fName);
+
+    // Assume non-collated (as fallback value).
+    // If everyone claims to have the file, use master to verify if
+    // collated is involved.
+
+    bool isCollated = false;
+
+    if (returnReduceAnd(found, UPstream::worldComm))
+    {
+        // Test for collated format.
+        // Note: can test only world-master. Since even host-collated will have
+        // same file format type for all processors
+        if (UPstream::master(UPstream::worldComm))
+        {
+            const bool oldParRun = UPstream::parRun(false);
+
+            if (IFstream is(fName); is.good())
+            {
+                IOobject io(name, instance, local, runTime);
+                io.readHeader(is);
+
+                isCollated = decomposedBlockData::isCollatedType(io);
+            }
+
+            UPstream::parRun(oldParRun);
+        }
+        Pstream::broadcast(isCollated, UPstream::worldComm);
+    }
+
+    // For collated, check that the corresponding blocks exist
+    if (isCollated)
+    {
+        found = checkFileExistenceCollated(handler, fName);
+    }
+
+
+    // Globally consistent information about who has the file
+    bitSet haveFileOnProc = bitSet::allGather(found, UPstream::worldComm);
+
+    if (verbose)
+    {
+        Info<< "Per processor availability of \""
+            << name << "\" file in " << instance/local << nl
+            << "    " << flatOutput(haveFileOnProc) << nl << endl;
+    }
+
+    return haveFileOnProc;
+}
+
+
+Foam::boolList Foam::haveMeshFile
+(
+    const word& name,           // eg "faces"
+    const fileName& instance,   // eg "constant"
+    const fileName& local,      // eg, polyMesh
+    const Time& runTime,
+    const bool verbose
+)
+{
+    const auto& handler = Foam::fileHandler();
+
+    const fileName fName
+    (
+        handler.filePath(runTime.path()/instance/local/name)
+    );
+
+    bool found = handler.isFile(fName);
+
+    // Assume non-collated (as fallback value).
+    // If everyone claims to have the file, use master to verify if
+    // collated is involved.
+
+    bool isCollated = false;
+
+    if (returnReduceAnd(found, UPstream::worldComm))
+    {
+        // Test for collated format.
+        // Note: can test only world-master. Since even host-collated will have
+        // same file format type for all processors
+        if (UPstream::master(UPstream::worldComm))
+        {
+            const bool oldParRun = UPstream::parRun(false);
+
+            if (IFstream is(fName); is.good())
+            {
+                IOobject io(name, instance, local, runTime);
+                io.readHeader(is);
+
+                isCollated = decomposedBlockData::isCollatedType(io);
+            }
+
+            UPstream::parRun(oldParRun);
+        }
+        Pstream::broadcast(isCollated, UPstream::worldComm);
+    }
+
+    // For collated, check that the corresponding blocks exist
+    if (isCollated)
+    {
+        found = checkFileExistenceCollated(handler, fName);
+    }
+
+
+    // Globally consistent information about who has a mesh
+    boolList haveFileOnProc
+    (
+        UPstream::allGatherValues<bool>(found, UPstream::worldComm)
+    );
+
+    if (verbose)
+    {
+        Info<< "Per processor availability of \""
+            << name << "\" file in " << instance/local << nl
+            << "    " << flatOutput(haveFileOnProc) << nl << endl;
+    }
+
+    return haveFileOnProc;
+}
+
 
 Foam::boolList Foam::haveMeshFile
 (
@@ -135,143 +302,49 @@ Foam::boolList Foam::haveMeshFile
     const bool verbose
 )
 {
-    #if 0
-
-    // Simple directory scanning - too fragile
-    bool found = checkFileExistence(runTime.path()/meshPath/meshFile);
-
-    #else
-
-    // Trimmed-down version of lookupAndCacheProcessorsPath
-    // with Foam::exists() check. No caching.
-
-    // Check for two conditions:
-    // - file has to exist
-    // - if collated the entry has to exist inside the file
-
-    // Note: bypass fileOperation::filePath(IOobject&) since has problems
-    //       with going to a different number of processors
-    //       (in collated format). Use file-based searching instead
-
     const auto& handler = Foam::fileHandler();
-    typedef fileOperation::procRangeType procRangeType;
 
     const fileName fName
     (
         handler.filePath(runTime.path()/meshPath/meshFile)
     );
+
     bool found = handler.isFile(fName);
-    if (returnReduceAnd(found)) // worldComm
+
+    // Assume non-collated (as fallback value).
+    // If everyone claims to have the file, use master to verify if
+    // collated is involved.
+
+    bool isCollated = false;
+
+    if (returnReduceAnd(found, UPstream::worldComm))
     {
-        // Bit tricky: avoid having all slaves open file since this involves
-        // reading it on master and broadcasting it. This fails if file > 2G.
-        // So instead only read on master
-
-        bool isCollated = false;
-
+        // Test for collated format.
         // Note: can test only world-master. Since even host-collated will have
         // same file format type for all processors
         if (UPstream::master(UPstream::worldComm))
         {
             const bool oldParRun = UPstream::parRun(false);
 
-            IFstream is(fName);
-            if (is.good())
+            if (IFstream is(fName); is.good())
             {
                 IOobject io(meshFile, meshPath, runTime);
                 io.readHeader(is);
 
                 isCollated = decomposedBlockData::isCollatedType(io);
             }
+
             UPstream::parRun(oldParRun);
         }
-        Pstream::broadcast(isCollated); //UPstream::worldComm
-
-
-        // Collect block-number in individual filenames (might differ
-        // on different processors)
-        if (isCollated)
-        {
-            const label nProcs = UPstream::nProcs(fileHandler().comm());
-            const label myProcNo = UPstream::myProcNo(fileHandler().comm());
-
-            // Collect file names on master of local communicator
-            const fileNameList fNames
-            (
-                Pstream::listGatherValues
-                (
-                    fName,
-                    fileHandler().comm(),
-                    UPstream::msgType()
-                )
-            );
-
-            // Collect local block number
-            label myBlockNumber = -1;
-            {
-                procRangeType group;
-                label proci = fileOperation::detectProcessorPath(fName, group);
-
-                if (proci == -1 && group.empty())
-                {
-                    // 'processorsXXX' format so contains all ranks
-                    // according to worldComm
-                    myBlockNumber = UPstream::myProcNo(UPstream::worldComm);
-                }
-                else
-                {
-                    // 'processorsXXX_n-m' format so check for the
-                    // relative rank
-                    myBlockNumber = myProcNo;
-                }
-            }
-            const labelList myBlockNumbers
-            (
-                Pstream::listGatherValues
-                (
-                    myBlockNumber,
-                    fileHandler().comm(),
-                    UPstream::msgType()
-                )
-            );
-
-
-
-            // Determine for all whether the filename exists in the collated
-            // file.
-            boolList allFound(nProcs, false);
-
-            if (UPstream::master(fileHandler().comm()))
-            {
-                // Store nBlocks and index of file that was used for nBlocks
-                label nBlocks = -1;
-                label blockRanki = -1;
-                forAll(fNames, ranki)
-                {
-                    if
-                    (
-                        blockRanki == -1
-                     || (fNames[ranki] != fNames[blockRanki])
-                    )
-                    {
-                        blockRanki = ranki;
-                        IFstream is(fNames[ranki]);
-                        nBlocks = decomposedBlockData::getNumBlocks(is);
-                    }
-
-                    allFound[ranki] = (myBlockNumbers[ranki] < nBlocks);
-                }
-            }
-
-            found = Pstream::listScatterValues
-            (
-                allFound,
-                fileHandler().comm(),
-                UPstream::msgType()
-            );
-        }
+        Pstream::broadcast(isCollated, UPstream::worldComm);
     }
-    #endif
+
+    // For collated, check that the corresponding blocks exist
+    if (isCollated)
+    {
+        found = checkFileExistenceCollated(handler, fName);
+    }
+
 
     // Globally consistent information about who has a mesh
     boolList haveFileOnProc

applications/utilities/parallelProcessing/redistributePar/loadOrCreateMesh.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012 OpenFOAM Foundation
-    Copyright (C) 2022-2023 OpenCFD Ltd.
+    Copyright (C) 2022-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -48,15 +48,32 @@ namespace Foam
 // Forward Declarations
 class faMesh;
 
-//- Check for availability of given file
-bool checkFileExistence(const fileName& fName);
+//- Check for availability of specified file (on all ranks)
+bitSet haveProcessorFile
+(
+    const word& name,           // eg, "faces"
+    const fileName& instance,   // eg, "constant"
+    const fileName& local,      // eg, "polyMesh"
+    const Time& runTime,
+    const bool verbose = true
+);
 
-//- Check for availability of specified mesh file (default: "faces")
+//- Check for availability of specified file (on all ranks)
+boolList haveMeshFile
+(
+    const word& name,           // eg, "faces"
+    const fileName& instance,   // eg, "constant"
+    const fileName& local,      // eg, "polyMesh"
+    const Time& runTime,
+    const bool verbose = true
+);
+
+//- Check for availability of specified mesh file
 boolList haveMeshFile
 (
     const Time& runTime,
     const fileName& meshPath,
-    const word& meshFile = "faces",
+    const word& meshFile,
     const bool verbose = true
 );
 
@@ -75,7 +92,7 @@ void masterMeshInstance
     fileName& pointsInstance
 );
 
-}
+} // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/utilities/postProcessing/dataConversion/foamToEnsight/checkFieldAvailability.H

@@ -21,44 +21,59 @@ Requires
 
 // Initially all possible objects that are available at the final time
 List<wordHashSet> availableRegionObjectNames(meshes.size());
-List<wordHashSet> availableFaRegionObjectNames(meshes.size());
+List<HashTable<wordHashSet>> availableFaRegionObjectNames(meshes.size());
 
 forAll(meshes, regioni)
 {
     const auto& mesh = meshes[regioni];
 
     IOobjectList objects;
-    IOobjectList faObjects;
+    HashTable<IOobjectList> faObjects;
 
     if (doConvertFields && !timeDirs.empty())
     {
+        // const word& checkTimeDir = timeDirs.back().name();
+
         // List of volume mesh objects for this instance
         objects = IOobjectList(mesh, timeDirs.back().name());
 
-        // List of area mesh objects (assuming single region)
-        faObjects = IOobjectList
-        (
-            mesh.time(),
-            timeDirs.back().name(),
-            faMesh::dbDir(mesh, word::null),
-            IOobjectOption::NO_REGISTER
-        );
-
         if (fieldSelector)
         {
             objects.filterObjects(fieldSelector);
-            faObjects.filterObjects(fieldSelector);
         }
 
         // Remove "*_0" restart fields
         objects.prune_0();
-        faObjects.prune_0();
 
         if (!doPointValues)
         {
             // Prune point fields if disabled
             objects.filterClasses(Foam::fieldTypes::is_point, true);
         }
+
+
+        // The finite-area regions and their fields for this volume region
+        // and instance
+        faObjects.reserve(areaRegionNames.size());
+        for (const word& areaName : areaRegionNames)
+        {
+            IOobjectList objs
+            (
+                faMesh::Registry(mesh),
+                timeDirs.back().name(),
+                polyMesh::regionName(areaName),
+                IOobjectOption::NO_REGISTER
+            );
+
+            if (fieldSelector)
+            {
+                objs.filterObjects(fieldSelector);
+            }
+
+            // Remove "*_0" restart fields
+            objs.prune_0();
+            faObjects(areaName) = std::move(objs);
+        }
     }
 
     // Volume fields
@@ -78,20 +93,29 @@ forAll(meshes, regioni)
     }
 
     // Area fields
-    if (!faObjects.empty())
+    for (const word& areaName : areaRegionNames)
     {
-        wordList objectNames(faObjects.sortedNames());
+        wordList objectNames;
 
-        // Check availability for all times... (assuming single region)
-        checkData
+        if (const auto iter = faObjects.cfind(areaName); iter.good())
+        {
+            objectNames = iter.val().sortedNames();
+
+            // Check availability for all times...
+            checkData
+            (
+                faMesh::Registry(mesh),
+                timeDirs,
+                objectNames,
+                polyMesh::regionName(areaName)
+            );
+        }
+
+        availableFaRegionObjectNames[regioni].emplace_set
         (
-            mesh.time(),
-            timeDirs,
-            objectNames,
-            faMesh::dbDir(mesh, word::null)
+            areaName,
+            std::move(objectNames)
         );
-
-        availableFaRegionObjectNames[regioni] = objectNames;
     }
 }
 

applications/utilities/postProcessing/dataConversion/foamToEnsight/createMeshAccounting.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021-2022 OpenCFD Ltd.
+    Copyright (C) 2021-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -15,12 +15,14 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
+// Cases and meshes per volume region
 PtrList<ensightCase> ensightCases(regionNames.size());
 PtrList<ensightMesh> ensightMeshes(regionNames.size());
 
-PtrList<faMesh> meshesFa(regionNames.size());
-PtrList<ensightCase> ensightCasesFa(regionNames.size());
-PtrList<ensightFaMesh> ensightMeshesFa(regionNames.size());
+// Per volume region can have multiple finite-area regions
+List<PtrDynList<faMesh>> meshesFa(regionNames.size());
+List<PtrDynList<ensightCase>> ensightCasesFa(regionNames.size());
+List<PtrDynList<ensightFaMesh>> ensightMeshesFa(regionNames.size());
 
 {
     forAll(regionNames, regioni)
@@ -45,6 +47,7 @@ PtrList<ensightFaMesh> ensightMeshesFa(regionNames.size());
             }
         }
 
+        // Volume mesh
         ensightMeshes.set
         (
             regioni,
@@ -59,33 +62,60 @@ PtrList<ensightFaMesh> ensightMeshesFa(regionNames.size());
             new ensightCase(ensCasePath, ensCaseName, caseOpts)
         );
 
-        if (doFiniteArea)
+        if (!doFiniteArea)
         {
-            autoPtr<faMesh> faMeshPtr(faMesh::TryNew(mesh));
+            continue;
+        }
+
+        // Note: not every volume region is guaranteed to have
+        // any or all area region(s)
+
+        meshesFa[regioni].reserve_exact(areaRegionNames.size());
+
+        for (const word& areaName : areaRegionNames)
+        {
+            auto faMeshPtr = faMesh::TryNew(areaName, mesh);
 
             if (faMeshPtr)
             {
-                ensightCasesFa.set
-                (
-                    regioni,
-                    new ensightCase
-                    (
-                        ensCasePath/"finite-area",
-                        "finite-area",
-                        caseOpts
-                    )
-                );
-
-                meshesFa.set(regioni, std::move(faMeshPtr));
-
-                ensightMeshesFa.set
-                (
-                    regioni,
-                    new ensightFaMesh(meshesFa[regioni])
-                );
-                ensightMeshesFa[regioni].verbose(optVerbose);
+                meshesFa[regioni].push_back(std::move(faMeshPtr));
             }
         }
+
+        // Setup the ensight components
+        const label nAreas = meshesFa[regioni].size();
+        ensightCasesFa[regioni].reserve_exact(nAreas);
+        ensightMeshesFa[regioni].reserve_exact(nAreas);
+
+        for (const faMesh& areaMesh : meshesFa[regioni])
+        {
+            // Use regionName() - automatically filters for defaultRegion
+            const word& areaName = areaMesh.regionName();
+
+            if (areaName.empty())
+            {
+                // single-region
+                ensightCasesFa[regioni].emplace_back
+                (
+                    ensCasePath/"finite-area",
+                    "finite-area",
+                    caseOpts
+                );
+            }
+            else
+            {
+                // multi-region
+                ensightCasesFa[regioni].emplace_back
+                (
+                    ensCasePath/"finite-area"/areaName,
+                    areaName,
+                    caseOpts
+                );
+            }
+
+            auto& ensFaMesh = ensightMeshesFa[regioni].emplace_back(areaMesh);
+            ensFaMesh.verbose(optVerbose);
+        }
     }
 }
 

applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C

@@ -124,6 +124,7 @@ Usage
 #include "OFstream.H"
 #include "Pstream.H"
 #include "HashOps.H"
+#include "PtrDynList.H"
 #include "regionProperties.H"
 
 #include "fvc.H"
@@ -172,6 +173,7 @@ int main(int argc, char *argv[])
     argList::addVerboseOption();
 
     #include "addAllRegionOptions.H"
+    #include "addAllFaRegionOptions.H"
 
     argList::addBoolOption
     (
@@ -439,6 +441,14 @@ int main(int argc, char *argv[])
     // Handle -allRegions, -regions, -region
     #include "getAllRegionOptions.H"
 
+    // Handle -all-area-regions, -area-regions, -area-region
+    #include "getAllFaRegionOptions.H"
+
+    if (!doFiniteArea)
+    {
+        areaRegionNames.clear();  // For consistency
+    }
+
     // ------------------------------------------------------------------------
     // Directory management
 
@@ -518,29 +528,32 @@ int main(int argc, char *argv[])
             polyMesh::readUpdateState meshState = mesh.readUpdate();
             const bool moving = (meshState != polyMesh::UNCHANGED);
 
-            // Ensight
+            // Ensight (fvMesh)
             auto& ensCase = ensightCases[regioni];
             auto& ensMesh = ensightMeshes[regioni];
 
-            // Finite-area (can be missing)
-            auto* ensFaCasePtr = ensightCasesFa.get(regioni);
-            auto* ensFaMeshPtr = ensightMeshesFa.get(regioni);
-
             ensCase.setTime(timeDirs[timei], timeIndex);
-            if (ensFaCasePtr)
+
+            // Finite-area (optional)
+            // Accounting exists for each volume region but may be empty
+            auto& ensFaCases = ensightCasesFa[regioni];
+            auto& ensFaMeshes = ensightMeshesFa[regioni];
+
+            for (auto& ensFaCase : ensFaCases)
             {
-                ensFaCasePtr->setTime(timeDirs[timei], timeIndex);
+                ensFaCase.setTime(timeDirs[timei], timeIndex);
             }
 
+            // Movement
             if (moving)
             {
                 ensMesh.expire();
                 ensMesh.correct();
 
-                if (ensFaMeshPtr)
+                for (auto& ensFaMesh : ensFaMeshes)
                 {
-                    ensFaMeshPtr->expire();
-                    ensFaMeshPtr->correct();
+                    ensFaMesh.expire();
+                    ensFaMesh.correct();
                 }
             }
 
@@ -555,15 +568,19 @@ int main(int argc, char *argv[])
                 }
 
                 // finite-area
-                if (ensFaCasePtr && ensFaMeshPtr)
+                forAll(ensFaMeshes, areai)
                 {
-                    autoPtr<ensightGeoFile> os =
-                        ensFaCasePtr->newGeometry(hasMovingMesh);
+                    const auto& ensFaCase = ensFaCases[areai];
+                    const auto& ensFaMesh = ensFaMeshes[areai];
 
-                    ensFaMeshPtr->write(os.ref());
+                    autoPtr<ensightGeoFile> os =
+                        ensFaCase.newGeometry(hasMovingMesh);
+
+                    ensFaMesh.write(os.ref());
                 }
             }
 
+
             // Objects at this time
             IOobjectList objects(mesh, runTime.timeName());
 
@@ -575,23 +592,37 @@ int main(int argc, char *argv[])
             // Volume, internal, point fields
             #include "convertVolumeFields.H"
 
-            // The finite-area objects at this time
-            IOobjectList faObjects;
-
-            if (ensFaMeshPtr)
+            // finite-area
+            forAll(ensFaMeshes, areai)
             {
-                faObjects =
-                    IOobjectList(ensFaMeshPtr->mesh(), runTime.timeName());
+                auto* ensFaCasePtr = ensFaCases.get(areai);
+                auto* ensFaMeshPtr = ensFaMeshes.get(areai);
 
-                faObjects.filterObjects
-                (
-                    availableFaRegionObjectNames[regioni]
-                );
+                // The finite-area region objects at this time
+                IOobjectList faObjects;
+
+                if (ensFaMeshPtr)
+                {
+                    const word& areaName = ensFaMeshPtr->mesh().name();
+
+                    faObjects = IOobjectList
+                    (
+                        faMesh::Registry(mesh),
+                        runTime.timeName(),
+                        polyMesh::regionName(areaName),
+                        IOobjectOption::NO_REGISTER
+                    );
+
+                    faObjects.filterObjects
+                    (
+                        availableFaRegionObjectNames[regioni](areaName)
+                    );
+                }
+
+                // The finiteArea fields
+                #include "convertAreaFields.H"
             }
 
-            // The finiteArea fields
-            #include "convertAreaFields.H"
-
             // Lagrangian fields
             #include "convertLagrangian.H"
         }
@@ -604,14 +635,13 @@ int main(int argc, char *argv[])
     // Write cases
     forAll(ensightCases, regioni)
     {
+        // finite-volume
         ensightCases[regioni].write();
-    }
 
-    forAll(ensightCasesFa, regioni)
-    {
-        if (ensightCasesFa.set(regioni))
+        // Finite-area (if any)
+        for (const auto& ensFaCase : ensightCasesFa[regioni])
         {
-            ensightCasesFa[regioni].write();
+            ensFaCase.write();
         }
     }
 

applications/utilities/postProcessing/dataConversion/foamToEnsight/readFields.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2023 OpenCFD Ltd.
+    Copyright (C) 2018-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -26,7 +26,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "readFields.H"
-#include "volFields.H"
 
 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 
@@ -42,26 +41,30 @@ Foam::label Foam::checkData
 
     wordHashSet goodFields;
 
+    IOobject io
+    (
+        "any-name",   // placeholder
+        "constant",   // placeholder
+        local,
+        obr,
+        IOobjectOption::NO_READ,
+        IOobjectOption::NO_WRITE,
+        IOobjectOption::NO_REGISTER
+    );
+
     for (const word& fieldName : objectNames)
     {
-        // // If prune_0() not previously used...
-        // if (objectNames.ends_with("_0")) continue;
+        // In case prune_0() not previously used...
+        if (fieldName.ends_with("_0")) continue;
 
         bool good = false;
 
         for (const instant& inst : timeDirs)
         {
-            good =
-                IOobject
-                (
-                    fieldName,
-                    inst.name(),
-                    local,
-                    obr,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE,
-                    IOobject::NO_REGISTER
-                ).typeHeaderOk<volScalarField>(false, false);
+            io.resetHeader(fieldName);
+            io.instance() = inst.name();
+
+            good = io.typeHeaderOk<regIOobject>(false, false);
 
             if (!good)
             {

applications/utilities/postProcessing/dataConversion/foamToVTK/convertAreaFields.H

@@ -22,24 +22,66 @@ Description
 //
 
 // No subsetting!
-if (doFiniteArea)
+if (doFiniteArea && !areaRegionNames.empty())
 {
     using reportFields = foamToVtkReportFields;
 
-    autoPtr<faMesh> faMeshPtr;
+    // The (region) polyMesh being used. No subsetting possible
+    const auto& basePolyMesh = meshProxy.baseMesh();
 
-    const label nAreaFields = faObjects.count(Foam::fieldTypes::is_area);
-
-    if (nAreaFields || withMeshIds)
+    for (const word& areaName : areaRegionNames)
     {
-        faMeshPtr = faMesh::TryNew(meshProxy.baseMesh());
-    }
+        const bool isDefaultRegion(polyMesh::regionName(areaName).empty());
 
-    if (faMeshPtr && (nAreaFields || withMeshIds))
-    {
-        const faMesh& areaMesh = faMeshPtr();
+        // CAUTION
+        // If we want to have constant access to the HashTable:
+        //
+        // (SEGFAULT)
+        // const auto& faObjs = faObjects.lookup(areaName, IOobjectList());
+        //
+        // Use an empty fallback to avoid binding to a temporary:
+        //
+        // const IOobjectList emptyObjectList;
+        // const auto& faObjs = faObjects.lookup(areaName, emptyObjectList);
 
-        reportFields::area(Info, faObjects);
+        // Since we do not need the area fields afterwards,
+        // just move them out from the HashTable
+
+        IOobjectList faObjs;
+        if (auto iter = faObjects.find(areaName); iter.good())
+        {
+            faObjs = std::move(iter.val());
+        }
+
+        const label nAreaFields = faObjs.count(Foam::fieldTypes::is_area);
+
+        autoPtr<faMesh> faMeshPtr;
+
+        if (nAreaFields || withMeshIds)
+        {
+            faMeshPtr = faMesh::TryNew(areaName, basePolyMesh);
+        }
+
+        if (!faMeshPtr)
+        {
+            if (!isDefaultRegion)
+            {
+                // Report any area region specified but missing
+                // - silently ignore region0
+                Info<< "No area-mesh [" << polyMesh::regionName(areaName)
+                    << "] on volume-region ["
+                    << basePolyMesh.regionName() << "]" << endl;
+            }
+            continue;
+        }
+        const auto& areaMesh = faMeshPtr();
+
+        Info<< "Using area-mesh [" << polyMesh::regionName(areaName)
+            << "] on volume-region ["
+            << basePolyMesh.regionName() << "]" << endl;
+
+
+        reportFields::area(Info, faObjs);
 
         const auto& pp = faMeshPtr->patch();
 
@@ -49,7 +91,12 @@ if (doFiniteArea)
             writeOpts,
             (
                 outputDir/regionDir/"finite-area"
-              / "finiteArea" + timeDesc
+              / (
+                   isDefaultRegion
+                 ? fileName("finiteArea")
+                 : fileName(areaName/areaName)
+                )
+              + timeDesc
             ),
             UPstream::parRun()
         );
@@ -96,7 +143,7 @@ if (doFiniteArea)
         (
             writer,
             areaMesh,
-            faObjects,
+            faObjs,
             true // syncPar
         );
 

applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK.C

@@ -135,13 +135,13 @@ Note
 #include "emptyPolyPatch.H"
 #include "volPointInterpolation.H"
 #include "faceZoneMesh.H"
+#include "faMesh.H"
 #include "areaFields.H"
 #include "fvMeshSubsetProxy.H"
 #include "faceSet.H"
 #include "pointSet.H"
 #include "HashOps.H"
 #include "regionProperties.H"
-#include "stringListOps.H"  // For stringListOps::findMatching()
 
 #include "Cloud.H"
 #include "readFields.H"
@@ -434,6 +434,7 @@ int main(int argc, char *argv[])
     argList::addOptionCompat("one-boundary", {"allPatches", 1806});
 
     #include "addAllRegionOptions.H"
+    #include "addAllFaRegionOptions.H"
 
     argList::addOption
     (
@@ -630,6 +631,16 @@ int main(int argc, char *argv[])
     // Handle -allRegions, -regions, -region
     #include "getAllRegionOptions.H"
 
+    // Handle -all-area-regions, -area-regions, -area-region
+    #include "getAllFaRegionOptions.H"
+
+    if (!doFiniteArea)
+    {
+        areaRegionNames.clear();  // For consistency
+    }
+
+    // ------------------------------------------------------------------------
+
     // Names for sets and zones
     word cellSelectionName;
     word faceSetName;
@@ -777,8 +788,11 @@ int main(int argc, char *argv[])
                 }
             }
 
+            // fvMesh fields
             IOobjectList objects;
-            IOobjectList faObjects;
+
+            // faMesh fields (multiple finite-area regions per volume region)
+            HashTable<IOobjectList> faObjects;
 
             if (doConvertFields)
             {
@@ -786,33 +800,51 @@ int main(int argc, char *argv[])
                 objects =
                     IOobjectList(meshProxy.baseMesh(), runTime.timeName());
 
-                // List of area mesh objects (assuming single region)
-                faObjects =
-                    IOobjectList
-                    (
-                        runTime,
-                        runTime.timeName(),
-                        faMesh::dbDir(meshProxy.baseMesh(), word::null),
-                        IOobjectOption::NO_REGISTER
-                    );
-
                 if (fieldSelector)
                 {
                     objects.filterObjects(fieldSelector);
-                    faObjects.filterObjects(fieldSelector);
                 }
-
                 // Remove "*_0" restart fields
                 objects.prune_0();
-                faObjects.prune_0();
 
                 if (!doPointValues)
                 {
                     // Prune point fields if disabled
                     objects.filterClasses(Foam::fieldTypes::is_point, true);
                 }
+
+
+                // Lists of finite-area fields
+                faObjects.reserve(areaRegionNames.size());
+
+                for (const word& areaName : areaRegionNames)
+                {
+                    // The finite-area objects for given area region.
+                    // Add as method to faMeshesRegistry?
+
+                    IOobjectList objs
+                    (
+                        faMesh::Registry(meshProxy.baseMesh()),
+                        runTime.timeName(),
+                        polyMesh::regionName(areaName),
+                        IOobjectOption::NO_REGISTER
+                    );
+
+                    if (fieldSelector)
+                    {
+                        objs.filterObjects(fieldSelector);
+                    }
+                    // Remove "*_0" restart fields
+                    objs.prune_0();
+
+                    if (!objs.empty())
+                    {
+                        faObjects.emplace_set(areaName, std::move(objs));
+                    }
+                }
             }
 
+
             if (processorFieldsOnly)
             {
                 // Processor-patches only and continue

applications/utilities/preProcessing/setFields/Make/options

@@ -1,10 +1,10 @@
 EXE_INC = \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/finiteArea/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/finiteArea/lnInclude
 
 EXE_LIBS = \
+    -lmeshTools \
     -lfiniteVolume \
     -lfiniteArea \
-    -lmeshTools \
     -lgenericPatchFields

applications/utilities/preProcessing/setFields/setFields.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2022-2024 OpenCFD Ltd.
+    Copyright (C) 2022-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -47,6 +47,7 @@ Description
 #include "areaFields.H"
 #include "coupledFvPatch.H"
 #include "coupledFaPatch.H"
+#include "regionProperties.H"
 
 using namespace Foam;
 
@@ -86,17 +87,18 @@ bool consumeUnusedType(const fieldDescription& fieldDesc, Istream& is)
     //? typedef GeometricField<Type, faePatchField, areaMesh> fieldType3;
     //? typedef GeometricField<Type, fvsPatchField, volMesh>  fieldType4;
 
-    if
+    const bool isExpectedType
     (
         fieldDesc.type() == fieldType1::typeName
      || fieldDesc.type() == fieldType2::typeName
-    )
+    );
+
+    if (isExpectedType)
     {
         (void) pTraits<Type>(is);
-        return true;
     }
 
-    return false;
+    return isExpectedType;
 }
 
 
@@ -122,13 +124,18 @@ bool setCellFieldType
 (
     const fieldDescription& fieldDesc,
     const fvMesh& mesh,
-    const labelList& selectedCells,
+    const labelUList& selectedCells,
     Istream& is
 )
 {
     typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
 
-    if (fieldDesc.type() != fieldType::typeName)
+    const bool isExpectedType
+    (
+        fieldDesc.type() == fieldType::typeName
+    );
+
+    if (!isExpectedType)
     {
         return false;
     }
@@ -143,7 +150,9 @@ bool setCellFieldType
         fieldDesc.name(),
         mesh.thisDb().time().timeName(),
         mesh.thisDb(),
-        IOobject::MUST_READ
+        IOobjectOption::MUST_READ,
+        IOobjectOption::NO_WRITE,
+        IOobjectOption::NO_REGISTER
     );
 
     bool found = fieldHeader.typeHeaderOk<fieldType>(true);
@@ -151,13 +160,8 @@ bool setCellFieldType
     if (!found)
     {
         // Fallback to "constant" directory
-        fieldHeader = IOobject
-        (
-            fieldDesc.name(),
-            mesh.thisDb().time().constant(),
-            mesh.thisDb(),
-            IOobject::MUST_READ
-        );
+        fieldHeader.resetHeader();
+        fieldHeader.instance() = mesh.thisDb().time().constant();
         found = fieldHeader.typeHeaderOk<fieldType>(true);
     }
 
@@ -171,7 +175,7 @@ bool setCellFieldType
 
         fieldType field(fieldHeader, mesh, false);
 
-        if (isNull(selectedCells) || selectedCells.size() == field.size())
+        if (isNull(selectedCells) || (selectedCells.size() == field.size()))
         {
             field.primitiveFieldRef() = fieldValue;
         }
@@ -204,7 +208,7 @@ bool setCellFieldType
             << "Field " << fieldDesc.name() << " not found" << endl;
     }
 
-    return true;
+    return isExpectedType;
 }
 
 
@@ -216,13 +220,18 @@ bool setAreaFieldType
 (
     const fieldDescription& fieldDesc,
     const faMesh& mesh,
-    const labelList& selectedFaces,
+    const labelUList& selectedFaces,
     Istream& is
 )
 {
     typedef GeometricField<Type, faPatchField, areaMesh> fieldType;
 
-    if (fieldDesc.type() != fieldType::typeName)
+    const bool isExpectedType
+    (
+        fieldDesc.type() == fieldType::typeName
+    );
+
+    if (!isExpectedType)
     {
         return false;
     }
@@ -237,7 +246,9 @@ bool setAreaFieldType
         fieldDesc.name(),
         mesh.thisDb().time().timeName(),
         mesh.thisDb(),
-        IOobject::MUST_READ
+        IOobjectOption::MUST_READ,
+        IOobjectOption::NO_WRITE,
+        IOobjectOption::NO_REGISTER
     );
 
     bool found = fieldHeader.typeHeaderOk<fieldType>(true);
@@ -245,13 +256,8 @@ bool setAreaFieldType
     if (!found)
     {
         // Fallback to "constant" directory
-        fieldHeader = IOobject
-        (
-            fieldDesc.name(),
-            mesh.thisDb().time().constant(),
-            mesh.thisDb(),
-            IOobject::MUST_READ
-        );
+        fieldHeader.resetHeader();
+        fieldHeader.instance() = mesh.thisDb().time().constant(),
         found = fieldHeader.typeHeaderOk<fieldType>(true);
     }
 
@@ -292,7 +298,7 @@ bool setAreaFieldType
             << "Field " << fieldDesc.name() << " not found" << endl;
     }
 
-    return true;
+    return isExpectedType;
 }
 
 
@@ -304,13 +310,18 @@ bool setFaceFieldType
 (
     const fieldDescription& fieldDesc,
     const fvMesh& mesh,
-    const labelList& selectedFaces,
+    const labelUList& selectedFaces,
     Istream& is
 )
 {
     typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
 
-    if (fieldDesc.type() != fieldType::typeName)
+    const bool isExpectedType
+    (
+        fieldDesc.type() == fieldType::typeName
+    );
+
+    if (!isExpectedType)
     {
         return false;
     }
@@ -325,7 +336,9 @@ bool setFaceFieldType
         fieldDesc.name(),
         mesh.thisDb().time().timeName(),
         mesh.thisDb(),
-        IOobject::MUST_READ
+        IOobjectOption::MUST_READ,
+        IOobjectOption::NO_WRITE,
+        IOobjectOption::NO_REGISTER
     );
 
     bool found = fieldHeader.typeHeaderOk<fieldType>(true);
@@ -333,13 +346,8 @@ bool setFaceFieldType
     if (!found)
     {
         // Fallback to "constant" directory
-        fieldHeader = IOobject
-        (
-            fieldDesc.name(),
-            mesh.thisDb().time().constant(),
-            mesh.thisDb(),
-            IOobject::MUST_READ
-        );
+        fieldHeader.resetHeader();
+        fieldHeader.instance() = mesh.thisDb().time().constant();
         found = fieldHeader.typeHeaderOk<fieldType>(true);
     }
 
@@ -355,15 +363,14 @@ bool setFaceFieldType
 
         // Create flat list of selected faces and their value.
         Field<Type> allBoundaryValues(mesh.nBoundaryFaces());
-        forAll(field.boundaryField(), patchi)
+        for (const auto& pfld : field.boundaryField())
         {
             SubField<Type>
             (
                 allBoundaryValues,
-                field.boundaryField()[patchi].size(),
-                field.boundaryField()[patchi].patch().start()
-              - mesh.nInternalFaces()
-            ) = field.boundaryField()[patchi];
+                pfld.size(),
+                pfld.patch().offset()
+            ) = pfld;
         }
 
         // Override
@@ -424,8 +431,7 @@ bool setFaceFieldType
                 (
                     allBoundaryValues,
                     fieldBf[patchi].size(),
-                    fieldBf[patchi].patch().start()
-                  - mesh.nInternalFaces()
+                    fieldBf[patchi].patch().offset()
                 );
             }
         }
@@ -445,7 +451,7 @@ bool setFaceFieldType
             << "Field " << fieldDesc.name() << " not found" << endl;
     }
 
-    return true;
+    return isExpectedType;
 }
 
 
@@ -460,7 +466,7 @@ struct setCellField
     (
         const fieldDescription& fieldDesc,
         const fvMesh& m,
-        const labelList& selectedCells,
+        const labelUList& selectedCells,
         Istream& is
     )
     {
@@ -477,11 +483,11 @@ struct setCellField
     class iNew
     {
         const fvMesh& mesh_;
-        const labelList& selected_;
+        const labelUList& selected_;
 
     public:
 
-        iNew(const fvMesh& mesh, const labelList& selectedCells)
+        iNew(const fvMesh& mesh, const labelUList& selectedCells) noexcept
         :
             mesh_(mesh),
             selected_(selectedCells)
@@ -528,7 +534,7 @@ struct setFaceField
     (
         const fieldDescription& fieldDesc,
         const fvMesh& m,
-        const labelList& selectedFaces,
+        const labelUList& selectedFaces,
         Istream& is
     )
     {
@@ -545,11 +551,11 @@ struct setFaceField
     class iNew
     {
         const fvMesh& mesh_;
-        const labelList& selected_;
+        const labelUList& selected_;
 
     public:
 
-        iNew(const fvMesh& mesh, const labelList& selectedFaces)
+        iNew(const fvMesh& mesh, const labelUList& selectedFaces) noexcept
         :
             mesh_(mesh),
             selected_(selectedFaces)
@@ -596,7 +602,7 @@ struct setAreaField
     (
         const fieldDescription& fieldDesc,
         const faMesh& m,
-        const labelList& selectedFaces,
+        const labelUList& selectedFaces,
         Istream& is
     )
     {
@@ -613,11 +619,11 @@ struct setAreaField
     class iNew
     {
         const faMesh& mesh_;
-        const labelList& selected_;
+        const labelUList& selected_;
 
     public:
 
-        iNew(const faMesh& mesh, const labelList& selectedFaces)
+        iNew(const faMesh& mesh, const labelUList& selectedFaces) noexcept
         :
             mesh_(mesh),
             selected_(selectedFaces)
@@ -675,6 +681,7 @@ int main(int argc, char *argv[])
     );
 
     #include "addRegionOption.H"
+    #include "addAllFaRegionOptions.H"
 
     // -------------------------
 
@@ -686,15 +693,59 @@ int main(int argc, char *argv[])
 
     #include "createNamedMesh.H"
 
-    autoPtr<faMesh> faMeshPtr;
+    // Handle -all-area-regions, -area-regions, -area-region
+    #include "getAllFaRegionOptions.H"
 
-    if (!args.found("no-finite-area"))
+    if (args.found("no-finite-area"))
     {
-        faMeshPtr = faMesh::TryNew(mesh);
+        areaRegionNames.clear();  // For consistency
     }
-    if (faMeshPtr)
+
+    // ------------------------------------------------------------------------
+
+    PtrList<faMesh> faMeshes;
+
+    // Setup all area meshes on this region
+    if (!areaRegionNames.empty())  // ie, !args.found("no-finite-area")
     {
-        Info<< "Detected finite-area mesh" << nl;
+        faMeshes.resize(areaRegionNames.size());
+
+        label nGoodRegions(0);
+
+        for (const word& areaName : areaRegionNames)
+        {
+            autoPtr<faMesh> faMeshPtr = faMesh::TryNew(areaName, mesh);
+
+            if (faMeshPtr)
+            {
+                faMeshes.set(nGoodRegions++, std::move(faMeshPtr));
+            }
+        }
+
+        faMeshes.resize(nGoodRegions);
+    }
+
+    if (faMeshes.size() == 1)
+    {
+        Info<< "Using finite-area mesh";
+        if
+        (
+            const word& name = polyMesh::regionName(faMeshes[0].name());
+            !name.empty()
+        )
+        {
+            Info<< " [" << name << "]";
+        }
+        Info<< nl;
+    }
+    else if (faMeshes.size() > 1)
+    {
+        Info<< "Detected finite-area meshes:";
+        for (const faMesh& areaMesh : faMeshes)
+        {
+            Info<< " [" << areaMesh.name() << "]";
+        }
+        Info<< nl;
     }
 
     const word dictName("setFieldsDict");
@@ -712,37 +763,45 @@ int main(int argc, char *argv[])
 
 
     // Default field values
+    if
+    (
+        const auto* eptr
+      = setFieldsDict.findEntry("defaultFieldValues", keyType::LITERAL)
+    )
     {
-        const entry* eptr =
-            setFieldsDict.findEntry("defaultFieldValues", keyType::LITERAL);
+        ITstream& is = eptr->stream();
 
-        if (eptr)
+        Info<< "Setting volume field default values" << endl;
+
+        PtrList<setCellField> defaultFieldValues
+        (
+            is,
+            setCellField::iNew(mesh, labelList::null())
+        );
+
+        for (const faMesh& areaMesh : faMeshes)
         {
-            ITstream& is = eptr->stream();
+            Info<< "Setting area field default values";
 
-            Info<< "Setting volume field default values" << endl;
-
-            PtrList<setCellField> defaultFieldValues
+            if
             (
-                is,
-                setCellField::iNew(mesh, labelList::null())
-            );
-
-            if (faMeshPtr)
+                const word& name = polyMesh::regionName(areaMesh.name());
+                !name.empty()
+            )
             {
-                const faMesh& areaMesh = faMeshPtr();
-                is.rewind();
-
-                Info<< "Setting area field default values" << endl;
-
-                PtrList<setAreaField> defaultFieldValues
-                (
-                    is,
-                    setAreaField::iNew(areaMesh, labelList::null())
-                );
+                Info<< " [" << name << "]";
             }
             Info<< endl;
+
+            is.rewind();
+
+            PtrList<setAreaField> defaultAreaFieldValues
+            (
+                is,
+                setAreaField::iNew(areaMesh, labelList::null())
+            );
         }
+        Info<< endl;
     }
 
 
@@ -768,11 +827,15 @@ int main(int argc, char *argv[])
 
             labelList selectedCells(subset.sortedToc());
 
-            Info<< "    Selected "
-                << returnReduce(selectedCells.size(), sumOp<label>())
-                << '/'
-                << returnReduce(mesh.nCells(), sumOp<label>())
-                << " cells" << nl;
+            {
+                FixedList<label, 2> stats;
+                stats[0] = selectedCells.size();
+                stats[1] = mesh.nCells();
+                reduce(stats, sumOp<label>());
+
+                Info<< "    Selected " << stats[0] << '/' << stats[1]
+                    << " cells" << nl;
+            }
 
             ITstream& is = region.dict().lookup("fieldValues");
 
@@ -806,10 +869,8 @@ int main(int argc, char *argv[])
                 setFaceField::iNew(mesh, selectedFaces)
             );
 
-            if (faMeshPtr)
+            for (const faMesh& areaMesh : faMeshes)
             {
-                const faMesh& areaMesh = faMeshPtr();
-
                 const labelUList& faceLabels = areaMesh.faceLabels();
 
                 // Transcribe from mesh faces to finite-area addressing
@@ -828,11 +889,15 @@ int main(int argc, char *argv[])
                 }
                 areaFaces.resize(nUsed);
 
-                Info<< "    Selected "
-                    << returnReduce(areaFaces.size(), sumOp<label>())
-                    << '/'
-                    << returnReduce(faceLabels.size(), sumOp<label>())
-                    << " area faces" << nl;
+                {
+                    FixedList<label, 2> stats;
+                    stats[0] = areaFaces.size();
+                    stats[1] = faceLabels.size();
+                    reduce(stats, sumOp<label>());
+
+                    Info<< "    Selected " << stats[0] << '/' << stats[1]
+                        << " area faces for " << areaMesh.name() << endl;
+                }
 
                 is.rewind();
 

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/faOptions/Make/files

@@ -9,7 +9,7 @@ faceSetOption/faceSetOption.C
 derivedSources=sources/derived
 
 $(derivedSources)/externalHeatFluxSource/externalHeatFluxSource.C
-$(derivedSources)/jouleHeatingSource/jouleHeatingSource.C
+$(derivedSources)/jouleHeatingSource/jouleHeatingSource.cxx
 $(derivedSources)/contactHeatFluxSource/contactHeatFluxSource.C
 $(derivedSources)/externalFileSource/externalFileSource.C
 

src/faOptions/faOption/faOption.C

@@ -81,6 +81,33 @@ Foam::fa::option::option
     active_(dict.getOrDefault("active", true)),
     log(true)
 {
+    // Suffix hint for variable names
+    if
+    (
+        coeffs_.readIfPresent("suffixing", suffixHint_)
+     || dict.readIfPresent("suffixing", suffixHint_)
+    )
+    {
+        Switch sw = Switch::find(suffixHint_);
+
+        if (sw.good())
+        {
+            if (!sw)  // No suffix
+            {
+                suffixHint_.clear();
+            }
+        }
+        else if (suffixHint_ == "default")
+        {
+            sw = true;
+        }
+
+        if (sw)  // Default suffix
+        {
+            suffixHint_ = '_' + regionName_;
+        }
+    }
+
     if (dict.readIfPresent("area", areaName_))
     {
         if (!sameRegionNames(areaName_, defaultAreaName))
@@ -95,9 +122,14 @@ Foam::fa::option::option
         }
     }
 
-    Log << incrIndent << indent << "Source: " << name_
-        << " [" << polyMesh::regionName(areaName_) << ']' << endl
-        << decrIndent;
+    Log << incrIndent << indent << "Source: " << name_;
+
+    if (!polyMesh::regionName(areaName_).empty())
+    {
+        Log<< " [" << areaName_ << ']';
+    }
+
+    Info<< decrIndent << endl;
 }
 
 
@@ -117,13 +149,20 @@ Foam::autoPtr<Foam::fa::option> Foam::fa::option::New
     coeffs.readIfPresent("area", areaName);
 
     Info<< indent
-        << "Selecting finite-area option, type " << modelType
-        << " [" << polyMesh::regionName(areaName) << ']';
+        << "Selecting finite-area option, type " << modelType;
 
-    if (!sameRegionNames(areaName, defaultAreaName))
+    if (sameRegionNames(areaName, defaultAreaName))
     {
-        Info<< " != " << defaultAreaName << nl;
+        if (!polyMesh::regionName(areaName).empty())
+        {
+            Info<< " [" << areaName << ']';
+        }
     }
+    else
+    {
+        Info<< " [" << areaName << "] != [" << defaultAreaName << ']';
+    }
+
     Info<< endl;
 
     mesh.time().libs().open

src/faOptions/faOption/faOption.H

@@ -65,10 +65,13 @@ Usage
                     --> the selected faOption settings | dictionary | no | -
       active   | Flag to (de)activate faOption         | bool  | no   | true
       log      | Flag to log faOption-related info     | bool  | no   | true
+      suffixing | Suffix hint for option variables     | word  | no   | -
     \endtable
 
 Note
   - Source/sink options are to be added to the right-hand side of equations.
+  - Suffixing (true|false|none|default|...) currently selects between
+    no suffix and ('_'+region).
 
 SourceFiles
     faOption.C
@@ -136,6 +139,9 @@ protected:
         //- The model region name (finite-area)
         word regionName_;
 
+        //- Suffix hint for variable names (default: "")
+        word suffixHint_;
+
 
     // Protected Member Functions
 
@@ -273,7 +279,7 @@ public:
         //- The source name
         const word& name() const noexcept { return name_; }
 
-        //- Return const access to the mesh database
+        //- Return const access to the volume mesh
         const fvMesh& mesh() const noexcept { return mesh_; }
 
         //- Return dictionary
@@ -304,6 +310,18 @@ public:
         inline bool active(bool on) noexcept;
 
 
+    // Helper Functions
+
+        //- The suffix hint when generating variable names
+        const word& suffixHint() const noexcept { return suffixHint_; }
+
+        //- Return the concatenation of \p base and the suffix hint
+        inline word suffixed(const char* base) const;
+
+        //- Return the concatenation of \p base and the suffix hint
+        inline word suffixed(const std::string& base) const;
+
+
     // Checks
 
         //- Is the source active?

src/faOptions/faOption/faOptionI.H

@@ -61,4 +61,18 @@ inline const Foam::volSurfaceMapping& Foam::fa::option::vsm() const
 }
 
 
+inline Foam::word
+Foam::fa::option::suffixed(const char* base) const
+{
+    return word(base + suffixHint_);
+}
+
+
+inline Foam::word
+Foam::fa::option::suffixed(const std::string& base) const
+{
+    return word(base + suffixHint_);
+}
+
+
 // ************************************************************************* //

src/faOptions/faOption/faOptions.C

@@ -27,7 +27,6 @@ License
 
 #include "faOptions.H"
 #include "faMesh.H"
-#include "faMeshesRegistry.H"
 #include "Time.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -68,7 +67,7 @@ Foam::IOobject createIOobject
         lookupName,
         mesh.time().constant(),
         // located under finite-area
-        faMeshesRegistry::New(mesh).thisDb(),
+        faMesh::Registry(mesh),
         IOobjectOption::MUST_READ,
         IOobjectOption::NO_WRITE,
         IOobjectOption::REGISTER
@@ -188,11 +187,7 @@ Foam::fa::options& Foam::fa::options::New
     );
 
     // Registered under finite-area?
-    auto* ptr =
-        faMeshesRegistry::New(mesh).thisDb().getObjectPtr<fa::options>
-        (
-            lookupName
-        );
+    auto* ptr = faMesh::Registry(mesh).getObjectPtr<fa::options>(lookupName);
 
     if (!ptr && polyMesh::regionName(defaultAreaName).empty())
     {

src/faOptions/sources/derived/contactHeatFluxSource/contactHeatFluxSource.C

@@ -30,7 +30,6 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "volFields.H"
 #include "famSup.H"
-#include "zeroGradientFaPatchFields.H"
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
@@ -57,13 +56,10 @@ Foam::fa::contactHeatFluxSource::contactHeatFluxSource
 :
     fa::faceSetOption(sourceName, modelType, dict, mesh, defaultAreaName),
     TName_(dict.getOrDefault<word>("T", "T")),
-    TprimaryName_(dict.get<word>("Tprimary")),
+    TprimaryName_(dict.getOrDefault<word>("Tprimary", "T")),
     Tprimary_(mesh_.lookupObject<volScalarField>(TprimaryName_)),
-    thicknessLayers_(),
-    kappaLayers_(),
     contactRes_(0),
-    curTimeIndex_(-1),
-    coupling_()
+    curTimeIndex_(-1)
 {
     fieldNames_.resize(1, TName_);
 
@@ -131,30 +127,30 @@ void Foam::fa::contactHeatFluxSource::addSup
     const label fieldi
 )
 {
-    if (isActive())
+    if (!isActive())
     {
-        DebugInfo
-            << name() << ": applying source to "
-            << eqn.psi().name() << endl;
+        return;
+    }
 
-        if (curTimeIndex_ != mesh().time().timeIndex())
-        {
-            tmp<DimensionedField<scalar, areaMesh>> htcw(htc());
+    DebugInfo<< name() << ": applying source to " << eqn.psi().name() << endl;
 
-            // Wall temperature - mapped from primary field to finite-area
-            auto Twall = DimensionedField<scalar, areaMesh>::New
-            (
-                "Tw_" + option::name(),
-                regionMesh(),
-                dimensionedScalar(dimTemperature, Zero)
-            );
+    if (curTimeIndex_ != mesh().time().timeIndex())
+    {
+        tmp<DimensionedField<scalar, areaMesh>> htcw(htc());
 
-            vsm().mapInternalToSurface<scalar>(Tprimary_, Twall.ref().field());
+        // Wall temperature - mapped from primary field to finite-area
+        auto Twall = DimensionedField<scalar, areaMesh>::New
+        (
+            "Tw_" + option::name(),
+            regionMesh(),
+            dimensionedScalar(dimTemperature, Zero)
+        );
 
-            eqn += -fam::Sp(htcw(), eqn.psi()) + htcw()*Twall;
+        vsm().mapInternalToSurface<scalar>(Tprimary_, Twall.ref().field());
 
-            curTimeIndex_ = mesh().time().timeIndex();
-        }
+        eqn += -fam::Sp(htcw(), eqn.psi()) + htcw()*Twall;
+
+        curTimeIndex_ = mesh().time().timeIndex();
     }
 }
 
@@ -190,8 +186,14 @@ bool Foam::fa::contactHeatFluxSource::read(const dictionary& dict)
 
         const labelList& patches = regionMesh().whichPolyPatches();
 
-        coupling_.clear();
-        coupling_.resize(patches.empty() ? 0 : (patches.last()+1));
+        if (patches.empty())
+        {
+            coupling_.clear();
+        }
+        else
+        {
+            coupling_.resize_null(patches.back()+1);
+        }
 
         for (const label patchi : patches)
         {

src/faOptions/sources/derived/contactHeatFluxSource/contactHeatFluxSource.H

@@ -40,6 +40,8 @@ Usage
     {
         // Mandatory entries (unmodifiable)
         type                contactHeatFluxSource;
+
+        // Optional entries (unmodifiable)
         Tprimary            <TprimaryFieldName>;
 
         // Optional entries (runtime modifiable)
@@ -60,12 +62,19 @@ Usage
     \table
       Property  | Description                       | Type  | Reqd | Dflt
       type      | Type name: contactHeatFluxSource  | word  | yes  | -
-      Tprimary  | Name of primary temperature field | word  | yes  | -
       T         | Name of operand temperature field | word  | no   | T
+      Tprimary  | Name of primary temperature field | word  | no   | T
       thicknessLayers | List of thicknesses of layers | scalarList | no | -
       kappaLayers | List of conductivities of layers | scalarList | cndtnl | -
     \endtable
 
+    Fields/variables used:
+    \table
+      Property  | Description               | Type   | Deflt
+      T         | Temperature               | shell  | T
+      Tprimary  | Temperature               | volume | T
+    \endtable
+
     The inherited entries are elaborated in:
      - \link faOption.H \endlink
      - \link faceSetOption.H \endlink
@@ -124,13 +133,13 @@ class contactHeatFluxSource
 
     // Private Data
 
-        //- Name of temperature field
+        //- Name of shell temperature field (default: "T")
         word TName_;
 
-        //- Name of primary temperature field
-        word TprimaryName_;
+        //- Name of volume temperature field (default: "T")
+        const word TprimaryName_;
 
-        //- Primary region temperature
+        //- Primary (volume) region temperature
         const volScalarField& Tprimary_;
 
         //- Thickness of layers

src/faOptions/sources/derived/externalHeatFluxSource/externalHeatFluxSource.C

@@ -72,10 +72,6 @@ Foam::fa::externalHeatFluxSource::externalHeatFluxSource
     fa::faceSetOption(sourceName, modelType, dict, m, defaultAreaName),
     mode_(operationModeNames.get("mode", dict)),
     TName_(dict.getOrDefault<word>("T", "T")),
-    Q_(nullptr),
-    q_(nullptr),
-    h_(nullptr),
-    Ta_(nullptr),
     emissivity_(dict.getOrDefault<scalar>("emissivity", 0))
 {
     fieldNames_.resize(1, TName_);
@@ -205,6 +201,11 @@ bool Foam::fa::externalHeatFluxSource::read(const dictionary& dict)
 
         mode_ = operationModeNames.get("mode", dict);
 
+        Q_.reset(nullptr);
+        q_.reset(nullptr);
+        h_.reset(nullptr);
+        Ta_.reset(nullptr);
+
         switch (mode_)
         {
             case fixedPower:

src/faOptions/sources/derived/externalHeatFluxSource/externalHeatFluxSource.H

@@ -103,7 +103,7 @@ Usage
     \verbatim
       power       | Use fixed power (supply Q)
       flux        | Use fixed heat flux (supply q)
-      coefficient | Use fixes heat transfer coefficient (supply h and T)
+      coefficient | Use fixes heat transfer coefficient (supply h and Ta)
     \endverbatim
 
 See also
@@ -159,7 +159,7 @@ private:
         //- Operation mode
         operationMode mode_;
 
-        //- Name of temperature field
+        //- Name of shell temperature field (default: "T")
         word TName_;
 
         //- Heat power [W]
@@ -174,7 +174,7 @@ private:
         //- Ambient temperature [K]
         autoPtr<Function1<scalar>> Ta_;
 
-        //- Optional surface emissivity for radiative transfer to ambient
+        //- Surface emissivity for radiative transfer to ambient (default: 0)
         scalar emissivity_;
 
 

src/faOptions/sources/derived/jouleHeatingSource/jouleHeatingSource.H

@@ -115,23 +115,30 @@ Usage
      - \link faceSetOption.H \endlink
 
 Note
-  - \c anisotropicElectricalConductivity=true enables
-  anisotropic (tensorial) electrical conductivity.
-  - \c anisotropicElectricalConductivity=false enables
-  isotropic (scalar) electrical conductivity.
-  - The electrical conductivity can be specified using either:
-    - If the \c sigma entry is present the electrical conductivity is specified
-    as a function of temperature using a \c Function1 type.
-    - If not present the \c sigma field will be read from file.
-    - If the \c anisotropicElectricalConductivity flag is set to \c true,
-    \c sigma should be specified as a tensor quantity.
+   If the \c sigma entry is present, the electrical conductivity is specified
+   as a function of temperature using a \c Function1 type, otherwise
+   the \c sigma field will be read from file.
+   When the \c anisotropicElectricalConductivity flag is set to \c true,
+   \c sigma should be specified as a \em tensor quantity instead of as
+   an isotropic \em scalar quantity.
+
+   BREAKING Naming changes from 2056 and earlier for the fields:
+   \table
+     Field   | Scoped names              | Scoped names (old)
+     V       | \<scope\>:V (suffix)      | \<scope\>:V_ + regionName
+     sigma   | \<scope\>:sigma (suffix)  | \<scope\>:sigma_ + regionName
+   \endtable
+
+   It is possible to replicate the older naming by specifying
+   the \c suffixing to ('_' + regionName).
 
 See also
   - Foam::Function1
+  - Foam::fv::jouleHeatingSource
 
 SourceFiles
-    jouleHeatingSource.C
-    jouleHeatingSourceTemplates.C
+    jouleHeatingSource.cxx
+    jouleHeatingSourceImpl.cxx
 
 \*---------------------------------------------------------------------------*/
 
@@ -188,13 +195,13 @@ class jouleHeatingSource
         void initialiseSigma
         (
             const dictionary& dict,
-            autoPtr<Function1<Type>>& sigmaVsTPtr
+            autoPtr<Function1<Type>>& sigmaFunctionPtr
         );
 
         //- Update the electrical conductivity field
         template<class Type>
         const GeometricField<Type, faPatchField, areaMesh>&
-        updateSigma(const autoPtr<Function1<Type>>& sigmaVsTPtr) const;
+        updateSigma(const autoPtr<Function1<Type>>& sigmaFunctionPtr) const;
 
 
 public:
@@ -229,22 +236,22 @@ public:
 
     // Member Functions
 
-        // Evaluation
+    // Evaluation
 
-            //- Add explicit contribution to compressible momentum equation
-            virtual void addSup
-            (
-                const areaScalarField& h,
-                const areaScalarField& rho,
-                faMatrix<scalar>& eqn,
-                const label fieldi
-            );
+        //- Add explicit contribution to energy equation
+        virtual void addSup
+        (
+            const areaScalarField& h,
+            const areaScalarField& rho,
+            faMatrix<scalar>& eqn,
+            const label fieldi
+        );
 
 
-        // IO
+    // IO
 
-            //- Read source dictionary
-            virtual bool read(const dictionary& dict);
+        //- Read source dictionary
+        virtual bool read(const dictionary& dict);
 };
 
 
@@ -255,12 +262,6 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#ifdef NoRepository
-    #include "jouleHeatingSourceTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //

src/faOptions/sources/derived/jouleHeatingSource/jouleHeatingSource.cxx

@@ -42,6 +42,12 @@ namespace fa
 }
 
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Implementation
+#include "jouleHeatingSourceImpl.cxx"
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::fa::jouleHeatingSource::jouleHeatingSource
@@ -59,7 +65,7 @@ Foam::fa::jouleHeatingSource::jouleHeatingSource
     (
         IOobject
         (
-            IOobject::scopedName(typeName, "V_" + regionName_),
+            suffixed(IOobject::scopedName(typeName, "V")),
             regionMesh().thisDb().time().timeName(),
             regionMesh().thisDb(),
             IOobject::MUST_READ,
@@ -68,8 +74,6 @@ Foam::fa::jouleHeatingSource::jouleHeatingSource
         ),
         regionMesh()
     ),
-    scalarSigmaVsTPtr_(nullptr),
-    tensorSigmaVsTPtr_(nullptr),
     curTimeIndex_(-1),
     nIter_(1),
     anisotropicElectricalConductivity_(false)
@@ -78,6 +82,8 @@ Foam::fa::jouleHeatingSource::jouleHeatingSource
 
     fa::option::resetApplied();
 
+    read(dict);
+
     if (anisotropicElectricalConductivity_)
     {
         Info<< "    Using tensor electrical conductivity" << endl;
@@ -90,8 +96,6 @@ Foam::fa::jouleHeatingSource::jouleHeatingSource
 
         initialiseSigma(coeffs_, scalarSigmaVsTPtr_);
     }
-
-    read(dict);
 }
 
 
@@ -105,12 +109,14 @@ void Foam::fa::jouleHeatingSource::addSup
     const label fieldi
 )
 {
-    if (isActive())
+    if (!isActive())
     {
-        DebugInfo
-            << name() << ": applying source to "
-            << eqn.psi().name() << endl;
+        return;
+    }
 
+    DebugInfo<< name() << ": applying source to " << eqn.psi().name() << endl;
+
+    {
         if (curTimeIndex_ != mesh().time().timeIndex())
         {
             for (label i = 0; i < nIter_; ++i)
@@ -118,8 +124,7 @@ void Foam::fa::jouleHeatingSource::addSup
                 if (anisotropicElectricalConductivity_)
                 {
                     // Update sigma as a function of T if required
-                    const areaTensorField& sigma =
-                        updateSigma(tensorSigmaVsTPtr_);
+                    const auto& sigma = updateSigma(tensorSigmaVsTPtr_);
 
                     // Solve the electrical potential equation
                     faScalarMatrix VEqn(fam::laplacian(h*sigma, V_));
@@ -129,8 +134,7 @@ void Foam::fa::jouleHeatingSource::addSup
                 else
                 {
                     // Update sigma as a function of T if required
-                    const areaScalarField& sigma =
-                        updateSigma(scalarSigmaVsTPtr_);
+                    const auto& sigma = updateSigma(scalarSigmaVsTPtr_);
 
                     // Solve the electrical potential equation
                     faScalarMatrix VEqn(fam::laplacian(h*sigma, V_));
@@ -145,7 +149,7 @@ void Foam::fa::jouleHeatingSource::addSup
         // Add the Joule heating contribution
         const word sigmaName
         (
-            IOobject::scopedName(typeName, "sigma_" + regionName_)
+            IOobject::scopedName(typeName, "sigma") + suffixHint()
         );
 
         areaVectorField gradV("gradV", fac::grad(V_));
@@ -162,15 +166,13 @@ void Foam::fa::jouleHeatingSource::addSup
 
         if (anisotropicElectricalConductivity_)
         {
-            const auto& sigma =
-                obr.lookupObject<areaTensorField>(sigmaName);
+            const auto& sigma = obr.lookupObject<areaTensorField>(sigmaName);
 
             tsource = (h*sigma & gradV) & gradV;
         }
         else
         {
-            const auto& sigma =
-                obr.lookupObject<areaScalarField>(sigmaName);
+            const auto& sigma = obr.lookupObject<areaScalarField>(sigmaName);
 
             tsource = (h*sigma*gradV) & gradV;
         }

src/faOptions/sources/derived/jouleHeatingSource/jouleHeatingSourceImpl.cxx

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2024 OpenCFD Ltd.
+    Copyright (C) 2019-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -33,16 +33,21 @@ template<class Type>
 void Foam::fa::jouleHeatingSource::initialiseSigma
 (
     const dictionary& dict,
-    autoPtr<Function1<Type>>& sigmaVsTPtr
+    autoPtr<Function1<Type>>& sigmaFunctionPtr
 )
 {
     typedef GeometricField<Type, faPatchField, areaMesh> FieldType;
 
+    const word sigmaName
+    (
+        IOobject::scopedName(typeName, "sigma") + suffixHint()
+    );
+
     auto& obr = regionMesh().thisDb();
 
     IOobject io
     (
-        IOobject::scopedName(typeName, "sigma_" + regionName_),
+        sigmaName,
         obr.time().timeName(),
         obr,
         IOobject::NO_READ,
@@ -52,11 +57,11 @@ void Foam::fa::jouleHeatingSource::initialiseSigma
 
     autoPtr<FieldType> tsigma;
 
-    if (dict.found("sigma"))
-    {
-        // Sigma to be defined using a Function1 type
-        sigmaVsTPtr = Function1<Type>::New("sigma", dict, &mesh_);
+    // Is sigma defined using a Function1 type?
+    sigmaFunctionPtr = Function1<Type>::NewIfPresent("sigma", dict, &mesh_);
 
+    if (sigmaFunctionPtr)
+    {
         tsigma.reset
         (
             new FieldType
@@ -89,31 +94,34 @@ template<class Type>
 const Foam::GeometricField<Type, Foam::faPatchField, Foam::areaMesh>&
 Foam::fa::jouleHeatingSource::updateSigma
 (
-    const autoPtr<Function1<Type>>& sigmaVsTPtr
+    const autoPtr<Function1<Type>>& sigmaFunctionPtr
 ) const
 {
     typedef GeometricField<Type, faPatchField, areaMesh> FieldType;
 
+    const word sigmaName
+    (
+        IOobject::scopedName(typeName, "sigma") + suffixHint()
+    );
+
     const auto& obr = regionMesh().thisDb();
 
-    auto& sigma =
-        obr.lookupObjectRef<FieldType>
-        (
-            IOobject::scopedName(typeName, "sigma_" + regionName_)
-        );
+    auto& sigma = obr.lookupObjectRef<FieldType>(sigmaName);
 
-    if (!sigmaVsTPtr)
+    if (!sigmaFunctionPtr)
     {
         // Electrical conductivity field, sigma, was specified by the user
         return sigma;
     }
 
+    const auto& sigmaFunction = sigmaFunctionPtr();
+
     const auto& T = obr.lookupObject<areaScalarField>(TName_);
 
     // Internal field
     forAll(sigma, i)
     {
-        sigma[i] = sigmaVsTPtr->value(T[i]);
+        sigma[i] = sigmaFunction.value(T[i]);
     }
 
 
@@ -127,7 +135,7 @@ Foam::fa::jouleHeatingSource::updateSigma
             const scalarField& Tbf = T.boundaryField()[patchi];
             forAll(pf, facei)
             {
-                pf[facei] = sigmaVsTPtr->value(Tbf[facei]);
+                pf[facei] = sigmaFunction.value(Tbf[facei]);
             }
         }
     }

src/finiteArea/faMesh/faMesh.C

@@ -28,6 +28,7 @@ License
 
 #include "faMesh.H"
 #include "faMeshBoundaryHalo.H"
+#include "faMeshesRegistry.H"
 #include "faGlobalMeshData.H"
 #include "Time.H"
 #include "polyMesh.H"
@@ -159,6 +160,12 @@ const Foam::objectRegistry* Foam::faMesh::registry(const polyMesh& pMesh)
 //     return obr.cfindObject<objectRegistry>(faMesh::prefix());
 // }
 
+// Forwarding
+const Foam::objectRegistry& Foam::faMesh::Registry(const polyMesh& pMesh)
+{
+    return faMeshesRegistry::Registry(pMesh);
+}
+
 
 const Foam::faMesh& Foam::faMesh::mesh
 (

src/finiteArea/faMesh/faMesh.H

@@ -751,10 +751,14 @@ public:
 
     // Database
 
-        //- The parent registry containing all finite-area meshes
-        //- on the polyMesh.
+        //- Find the singleton parent registry (on the polyMesh)
+        //- that contains all objects related to finite-area.
         static const objectRegistry* registry(const polyMesh& pMesh);
 
+        //- Return the singleton parent registry (on the polyMesh)
+        //- that contains all objects related to finite-area.
+        static const objectRegistry& Registry(const polyMesh& pMesh);
+
         //- The single-region finite-area region on the polyMesh.
         //- Uses lookupObject semantics - Fatal if non-existent
         static const faMesh& mesh(const polyMesh& pMesh);

src/finiteArea/faMesh/faMeshRegistry.C

@@ -27,8 +27,8 @@ License
 
 #include "faMesh.H"
 #include "faMeshesRegistry.H"
-#include "Time.H"
 #include "polyMesh.H"
+#include "Time.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -43,6 +43,7 @@ Foam::faMeshRegistry::faMeshRegistry
         IOobject
         (
             (areaName.empty() ? polyMesh::defaultRegion : areaName),
+            mesh.thisDb().time().timeName(),
             faMeshesRegistry::New(mesh).thisDb(),
             IOobjectOption::NO_READ,
             IOobjectOption::AUTO_WRITE,

src/finiteArea/faMesh/faMeshesRegistry.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2023-2024 OpenCFD Ltd.
+    Copyright (C) 2023-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -86,6 +86,8 @@ bool Foam::faMeshesRegistry::writeObject
     const bool writeOnProc
 ) const
 {
+    // Could also restrict to faMesh only...
+    //
     // for (const faMesh& m : objects_.csorted<faMesh>())
     // {
     //     m.write(writeOnProc);

src/finiteArea/faMesh/faMeshesRegistry.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2023-2024 OpenCFD Ltd.
+    Copyright (C) 2023-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -144,14 +144,34 @@ public:
         explicit faMeshesRegistry(const polyMesh& mesh);
 
 
+    // Factory Methods
+
+        //- Return the registry of objects on the singleton.
+        //  Same as New(mesh).thisDb()
+        FOAM_NO_DANGLING_REFERENCE  //< Reference stored in registry
+        static const objectRegistry& Registry(const polyMesh& mesh)
+        {
+            return MeshObject_type::New(mesh).thisDb();
+        }
+
+
     // Database
 
-        //- Return the object registry
+        // It redirects to the private objects but uses some
+        // objectRegistry method naming
+
+        //- The registry of the objects
         const objectRegistry& thisDb() const noexcept
         {
             return objects_;
         }
 
+        //- Local relative to time
+        const fileName& dbDir() const
+        {
+            return objects_.dbDir();
+        }
+
         //- The polyMesh reference
         const polyMesh& mesh() const noexcept
         {

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/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/fvOptions/Make/files

@@ -14,7 +14,7 @@ $(derivedSources)/buoyancyEnergy/buoyancyEnergy.C
 $(derivedSources)/buoyancyForce/buoyancyForce.C
 $(derivedSources)/directionalPressureGradientExplicitSource/directionalPressureGradientExplicitSource.C
 $(derivedSources)/explicitPorositySource/explicitPorositySource.C
-$(derivedSources)/jouleHeatingSource/jouleHeatingSource.C
+$(derivedSources)/jouleHeatingSource/jouleHeatingSource.cxx
 $(derivedSources)/meanVelocityForce/meanVelocityForce.C
 $(derivedSources)/meanVelocityForce/patchMeanVelocityForce/patchMeanVelocityForce.C
 $(derivedSources)/multiphaseStabilizedTurbulence/multiphaseStabilizedTurbulence.C

src/fvOptions/sources/derived/jouleHeatingSource/jouleHeatingSource.H

@@ -137,25 +137,28 @@ Note
   anisotropic (vectorial) electrical conductivity.
   - \c anisotropicElectricalConductivity=false enables
   isotropic (scalar) electrical conductivity.
-  - The electrical conductivity can be specified using either:
+
+  The electrical conductivity can be specified using either:
     - If the \c sigma entry is present the electrical conductivity is specified
     as a function of temperature using a \c Function1 type.
     - If not present the \c sigma field will be read from file.
     - If the \c anisotropicElectricalConductivity flag is set to \c true,
     \c sigma should be specified as a vector quantity.
+  .
 
 See also
   - Foam::Function1
   - Foam::coordSystem
+  - Foam::fa::jouleHeatingSource
 
 SourceFiles
-    jouleHeatingSource.C
-    jouleHeatingSourceTemplates.C
+    jouleHeatingSource.cxx
+    jouleHeatingSourceImpl.cxx
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef fv_jouleHeatingSource_H
-#define fv_jouleHeatingSource_H
+#ifndef Foam_fv_jouleHeatingSource_H
+#define Foam_fv_jouleHeatingSource_H
 
 #include "fvOption.H"
 #include "Function1.H"
@@ -185,9 +188,6 @@ class jouleHeatingSource
         //- Electrical potential field / [V]
         volScalarField V_;
 
-        //- Flag to indicate that the electrical conductivity is anisotropic
-        bool anisotropicElectricalConductivity_;
-
         //- Electrical conductivity as a scalar function of temperature
         autoPtr<Function1<scalar>> scalarSigmaVsTPtr_;
 
@@ -200,6 +200,9 @@ class jouleHeatingSource
         //- Current time index (used for updating)
         label curTimeIndex_;
 
+        //- Flag to indicate that the electrical conductivity is anisotropic
+        bool anisotropicElectricalConductivity_;
+
 
     // Private Member Functions
 
@@ -214,13 +217,13 @@ class jouleHeatingSource
         void initialiseSigma
         (
             const dictionary& dict,
-            autoPtr<Function1<Type>>& sigmaVsTPtr
+            autoPtr<Function1<Type>>& sigmaFunctionPtr
         );
 
         //- Update the electrical conductivity field
         template<class Type>
         const GeometricField<Type, fvPatchField, volMesh>&
-        updateSigma(const autoPtr<Function1<Type>>& sigmaVsTPtr) const;
+        updateSigma(const autoPtr<Function1<Type>>& sigmaFunctionPtr) const;
 
 
 public:
@@ -253,21 +256,21 @@ public:
 
     // Member Functions
 
-        // Evaluation
+    // Evaluation
 
-            //- Add explicit contribution to compressible momentum equation
-            virtual void addSup
-            (
-                const volScalarField& rho,
-                fvMatrix<scalar>& eqn,
-                const label fieldi
-            );
+        //- Add explicit contribution to energy equation
+        virtual void addSup
+        (
+            const volScalarField& rho,
+            fvMatrix<scalar>& eqn,
+            const label fieldi
+        );
 
 
-        // IO
+    // IO
 
-            //- Read source dictionary
-            virtual bool read(const dictionary& dict);
+        //- Read source dictionary
+        virtual bool read(const dictionary& dict);
 };
 
 
@@ -278,12 +281,6 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#ifdef NoRepository
-    #include "jouleHeatingSourceTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //

src/fvOptions/sources/derived/jouleHeatingSource/jouleHeatingSource.cxx

@@ -29,7 +29,6 @@ License
 #include "fvMatrices.H"
 #include "fvmLaplacian.H"
 #include "fvcGrad.H"
-#include "zeroGradientFvPatchField.H"
 #include "basicThermo.H"
 #include "addToRunTimeSelectionTable.H"
 
@@ -45,6 +44,12 @@ namespace fv
 }
 
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Implementation
+#include "jouleHeatingSourceImpl.cxx"
+
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 Foam::tmp<Foam::volSymmTensorField>
@@ -115,11 +120,8 @@ Foam::fv::jouleHeatingSource::jouleHeatingSource
         ),
         mesh
     ),
-    anisotropicElectricalConductivity_(false),
-    scalarSigmaVsTPtr_(nullptr),
-    vectorSigmaVsTPtr_(nullptr),
-    csysPtr_(nullptr),
-    curTimeIndex_(-1)
+    curTimeIndex_(-1),
+    anisotropicElectricalConductivity_(false)
 {
     // Set the field name to that of the energy
     // field from which the temperature is obtained
@@ -162,7 +164,7 @@ void Foam::fv::jouleHeatingSource::addSup
         else
         {
             // Update sigma as a function of T if required
-            const volScalarField& sigma = updateSigma(scalarSigmaVsTPtr_);
+            const auto& sigma = updateSigma(scalarSigmaVsTPtr_);
 
             // Solve the electrical potential equation
             fvScalarMatrix VEqn(fvm::laplacian(sigma, V_));
@@ -226,7 +228,7 @@ bool Foam::fv::jouleHeatingSource::read(const dictionary& dict)
 
             initialiseSigma(coeffs_, scalarSigmaVsTPtr_);
 
-            csysPtr_.clear();  // Do not need coordinate system
+            csysPtr_.reset(nullptr);  // Do not need coordinate system
         }
 
         return true;

src/fvOptions/sources/derived/jouleHeatingSource/jouleHeatingSourceImpl.cxx

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2024 OpenCFD Ltd.
+    Copyright (C) 2016-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,18 +27,25 @@ License
 
 #include "emptyFvPatchField.H"
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 template<class Type>
 void Foam::fv::jouleHeatingSource::initialiseSigma
 (
     const dictionary& dict,
-    autoPtr<Function1<Type>>& sigmaVsTPtr
+    autoPtr<Function1<Type>>& sigmaFunctionPtr
 )
 {
     typedef GeometricField<Type, fvPatchField, volMesh> FieldType;
 
+    const word sigmaName
+    (
+        IOobject::scopedName(typeName, "sigma")
+    );
+
     IOobject io
     (
-        IOobject::scopedName(typeName, "sigma"),
+        sigmaName,
         mesh_.time().timeName(),
         mesh_.thisDb(),
         IOobject::NO_READ,
@@ -48,11 +55,11 @@ void Foam::fv::jouleHeatingSource::initialiseSigma
 
     autoPtr<FieldType> tsigma;
 
-    if (dict.found("sigma"))
-    {
-        // Sigma to be defined using a Function1 type
-        sigmaVsTPtr = Function1<Type>::New("sigma", dict, &mesh_);
+    // Is sigma defined using a Function1 ?
+    sigmaFunctionPtr = Function1<Type>::NewIfPresent("sigma", dict, &mesh_);
 
+    if (sigmaFunctionPtr)
+    {
         tsigma.reset
         (
             new FieldType
@@ -85,28 +92,32 @@ template<class Type>
 const Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>&
 Foam::fv::jouleHeatingSource::updateSigma
 (
-    const autoPtr<Function1<Type>>& sigmaVsTPtr
+    const autoPtr<Function1<Type>>& sigmaFunctionPtr
 ) const
 {
     typedef GeometricField<Type, fvPatchField, volMesh> FieldType;
 
-    auto& sigma = mesh_.lookupObjectRef<FieldType>
+    const word sigmaName
     (
         IOobject::scopedName(typeName, "sigma")
     );
 
-    if (!sigmaVsTPtr)
+    auto& sigma = mesh_.lookupObjectRef<FieldType>(sigmaName);
+
+    if (!sigmaFunctionPtr)
     {
         // Electrical conductivity field, sigma, was specified by the user
         return sigma;
     }
 
+    const auto& sigmaFunction = sigmaFunctionPtr();
+
     const auto& T = mesh_.lookupObject<volScalarField>(TName_);
 
     // Internal field
     forAll(sigma, i)
     {
-        sigma[i] = sigmaVsTPtr->value(T[i]);
+        sigma[i] = sigmaFunction.value(T[i]);
     }
 
 
@@ -120,7 +131,7 @@ Foam::fv::jouleHeatingSource::updateSigma
             const scalarField& Tbf = T.boundaryField()[patchi];
             forAll(pf, facei)
             {
-                pf[facei] = sigmaVsTPtr->value(Tbf[facei]);
+                pf[facei] = sigmaFunction.value(Tbf[facei]);
             }
         }
     }

src/parallel/reconstruct/faReconstruct/faMeshReconstructor.cxx

@@ -480,6 +480,7 @@ void Foam::faMeshReconstructor::createMesh()
     (
         new faMesh
         (
+            procMesh_.name(),
             *serialVolMesh_,
             identity(singlePatchFaces_.size())  // faceLabels
         )

src/regionFaModels/KirchhoffShell/KirchhoffShell.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2023 OpenCFD Ltd.
+    Copyright (C) 2019-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -144,11 +144,23 @@ KirchhoffShell::KirchhoffShell
 )
 :
     vibrationShellModel(modelType, mesh, dict),
+    h_
+    (
+        IOobject
+        (
+            dict.getOrDefault<word>("h", suffixed("hs")),
+            regionMesh().time().timeName(),
+            regionMesh().thisDb(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        regionMesh()
+    ),
     ps_
     (
         IOobject
         (
-            "ps_" + regionName_,
+            dict.getOrDefault<word>("ps", suffixed("ps")),
             regionMesh().time().timeName(),
             regionMesh().thisDb(),
             IOobject::READ_IF_PRESENT,
@@ -157,23 +169,11 @@ KirchhoffShell::KirchhoffShell
         regionMesh(),
         dimensionedScalar(dimPressure, Zero)
     ),
-    h_
-    (
-        IOobject
-        (
-            "h_" + regionName_,
-            regionMesh().time().timeName(),
-            regionMesh().thisDb(),
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        regionMesh()
-    ),
     laplaceW_
     (
         IOobject
         (
-            "laplaceW_" + regionName_,
+            suffixed("laplaceW"),
             regionMesh().time().timeName(),
             regionMesh().thisDb(),
             IOobject::NO_READ,
@@ -186,7 +186,7 @@ KirchhoffShell::KirchhoffShell
     (
         IOobject
         (
-            "laplace2W_" + regionName_,
+            suffixed("laplace2W"),
             regionMesh().time().timeName(),
             regionMesh().thisDb(),
             IOobject::NO_READ,
@@ -199,7 +199,7 @@ KirchhoffShell::KirchhoffShell
     (
         IOobject
         (
-            "w0_" + regionName_,
+            suffixed("w0"),
             regionMesh().time().timeName(),
             regionMesh().thisDb(),
             IOobject::NO_READ,
@@ -212,7 +212,7 @@ KirchhoffShell::KirchhoffShell
     (
         IOobject
         (
-            "w00_" + regionName_,
+            suffixed("w00"),
             regionMesh().time().timeName(),
             regionMesh().thisDb(),
             IOobject::NO_READ,
@@ -225,7 +225,7 @@ KirchhoffShell::KirchhoffShell
     (
         IOobject
         (
-            "laplaceW0_" + regionName_,
+            suffixed("laplaceW0"),
             regionMesh().time().timeName(),
             regionMesh().thisDb(),
             IOobject::NO_READ,
@@ -238,7 +238,7 @@ KirchhoffShell::KirchhoffShell
     (
         IOobject
         (
-            "laplace2W0_" + regionName_,
+            suffixed("laplace2W0"),
             regionMesh().time().timeName(),
             regionMesh().thisDb(),
             IOobject::NO_READ,
@@ -256,6 +256,7 @@ KirchhoffShell::KirchhoffShell
     init(dict);
 }
 
+
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 void KirchhoffShell::preEvolveRegion()

src/regionFaModels/KirchhoffShell/KirchhoffShell.H

@@ -24,7 +24,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::regionFaModels::KirchhoffShell
+    Foam::regionModels::KirchhoffShell (Foam::regionFaModels)
 
 Description
     Vibration-shell finite-area model.
@@ -51,9 +51,25 @@ Usage
     \table
       Property   | Description                         | Type  | Reqd | Deflt
       vibrationShellModel | Type name: KirchhoffShell  | word  | yes  | -
-      f0         | Damping coefficient [1/s]           | scalar | yes | -
-      f1         | Damping coefficient [1/s]           | scalar | yes | -
-      f2         | Damping coefficient [1/s]           | scalar | yes | -
+      f0         | Damping coefficient [1/s]       | scalar | yes | -
+      f1         | Damping coefficient [1/s]       | scalar | yes | -
+      f2         | Damping coefficient [1/s]       | scalar | yes | -
+      h          | Name of thickness field         | word  | no   | hs (suffix)
+      ps         | Name of pressure field          | word  | no   | ps (suffix)
+    \endtable
+
+    Fields/variables used:
+    \table
+      Property  | Description               | Type  | Deflt
+      h         | Thickness                 | shell | hs (suffix)
+      ps        | Pressure on shell         | shell | ps (suffix)
+    \endtable
+
+    Naming changes from 2056 and earlier:
+    \table
+      Keyword  | Description               | Keyword (old) | Deflt (old)
+      h        | Thickness (shell)         | -   | "h_" + regionName
+      ps       | Pressure (shell)          | -   | "ps_" + regionName
     \endtable
 
     The inherited entries are elaborated in:
@@ -90,12 +106,12 @@ class KirchhoffShell
 
         // Source term fields
 
+        //- Shell thickness [m]
+        areaScalarField h_;
+
         //- External surface source [Pa]
         const areaScalarField ps_;
 
-        //- Thickness [m]
-        areaScalarField h_;
-
         //- Laplace of the displacement
         areaScalarField laplaceW_;
 

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

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

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

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

src/regionFaModels/regionFaModel/regionFaModel.C

@@ -27,7 +27,6 @@ License
 
 #include "regionFaModel.H"
 #include "faMesh.H"
-#include "faMeshesRegistry.H"
 #include "Time.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -74,7 +73,7 @@ Foam::IOobject createModelIOobject
     (
         objName,
         mesh.time().constant(),
-        faMeshesRegistry::New(mesh).thisDb(),
+        faMesh::Registry(mesh),
         IOobjectOption::NO_READ,
         IOobjectOption::NO_WRITE,
         IOobjectOption::REGISTER
@@ -141,12 +140,71 @@ Foam::IOobject createPropertiesIOobject
 
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
-void Foam::regionModels::regionFaModel::constructMeshObjects()
+void Foam::regionModels::regionFaModel::constructMeshObjects
+(
+    // Just for error reference
+    const dictionary& dict
+)
 {
+    regionMeshPtr_.reset(nullptr);
+
+    #if 1
     regionMeshPtr_.reset
     (
         new faMesh(areaName_, primaryMesh_)
     );
+
+    #else
+
+    // With try/catch and error messages
+
+    // DIY
+    // regionMeshPtr_ = faMesh::TryNew(areaName_, primaryMesh_);
+
+    // More heavy handed, but gives a better chance of locating
+    // the source of the error.
+    {
+        const bool oldThrowingError = FatalError.throwing(true);
+        const bool oldThrowingIOerr = FatalIOError.throwing(true);
+
+        try
+        {
+            regionMeshPtr_.reset
+            (
+                new faMesh(areaName_, primaryMesh_)
+            );
+        }
+        catch (const Foam::error& err)
+        {
+            Warning << err << nl << endl;
+
+            // Trickery to get original message
+            err.write(Warning, false);
+        }
+        catch (const Foam::IOerror& err)
+        {
+            Warning << err << nl << endl;
+
+            // Trickery to get original message
+            err.write(Warning, false);
+        }
+
+        FatalError.throwing(oldThrowingError);
+        FatalIOError.throwing(oldThrowingIOerr);
+    }
+
+    if (!regionMeshPtr_)
+    {
+        FatalError
+            << "Failed to create finite-area mesh [" << areaName_
+            << "] for model: "<< modelName_ << nl
+            << "A common cause is an incorrect or"
+               " missing 'area' entry in the setup" << nl
+            << ">>>>" << nl
+            << dict.relativeName() << dict << "<<<<" << endl
+            << exit(FatalError);
+    }
+    #endif
 }
 
 
@@ -234,7 +292,35 @@ Foam::regionModels::regionFaModel::regionFaModel
     regionName_(dict.get<word>("region")),
     coeffs_(dict.subOrEmptyDict(modelName + "Coeffs"))
 {
-    constructMeshObjects();
+    // Suffix hint for variable names
+    if
+    (
+        coeffs_.readIfPresent("suffixing", suffixHint_)
+     || dict.readIfPresent("suffixing", suffixHint_)
+    )
+    {
+        Switch sw = Switch::find(suffixHint_);
+
+        if (sw.good())
+        {
+            if (!sw)  // No suffix
+            {
+                suffixHint_.clear();
+            }
+        }
+        else if (suffixHint_ == "default")
+        {
+            // Use default suffix
+            sw = true;
+        }
+
+        if (sw)  // Default suffix
+        {
+            suffixHint_ = '_' + regionName_;
+        }
+    }
+
+    constructMeshObjects(dict);
     initialise();
 
     if (readFields)
@@ -250,9 +336,14 @@ void Foam::regionModels::regionFaModel::evolve()
 {
     if (active_)
     {
-        Info<< "\nEvolving " << modelName_ << " for region "
-            << regionMesh().name() << " : "
-            << polyMesh::regionName(areaName_) << endl;
+        Info<< "\nEvolving " << modelName_
+            << " for region " << regionMesh().name();
+
+        if (!polyMesh::regionName(areaName_).empty())
+        {
+            Info<< " [" << areaName_ << "]";
+        }
+        Info<< endl;
 
         preEvolveRegion();
 
@@ -265,7 +356,7 @@ void Foam::regionModels::regionFaModel::evolve()
         {
             Info<< incrIndent;
             info();
-            Info<< endl << decrIndent;
+            Info<< decrIndent << endl;
         }
     }
 }

src/regionFaModels/regionFaModel/regionFaModel.H

@@ -57,9 +57,20 @@ Usage
       region     | Name of operand region              | word  | yes  | -
       area       | Name of the finite-area mesh        | word  | no   | region0
       active     | Flag to activate the model          | bool  | yes  | -
+      suffixing  | Suffix hint for model variables     | word  | no   | -
       infoOutput | Flag to activate information output | bool  | no   | false
     \endtable
 
+Note
+    The \c suffixing parameter is a hint that individual models may use
+    when automatically generating variable names internally. For example,
+    a model \em may provide ("T" + suffixing) and ("q" + suffixing)
+    as its default names temperature and flux fields, respectively.
+    If the user specifies \c suffixing = "_foo", those \em default names
+    would then become "T_foo" and "q_foo", respectively.
+    Suffixing (true|false|none|default|...) currently selects between
+    no suffix and ('_'+region).
+
 SourceFiles
     regionFaModelI.H
     regionFaModel.C
@@ -90,7 +101,7 @@ class regionFaModel
     // Private Member Functions
 
         //- Construct region mesh and fields
-        void constructMeshObjects();
+        void constructMeshObjects(const dictionary&);
 
         //- Initialise the region
         void initialise();
@@ -118,18 +129,21 @@ protected:
         //- Model name
         const word modelName_;
 
-        //- The finite-area mesh name
+        //- Suffix hint for automatic model variable names (default: "")
+        word suffixHint_;
+
+        //- The finite-area mesh name (default: region0)
         word areaName_;
 
         //- Region name
         word regionName_;
 
-        //- Pointer to the region mesh database
-        autoPtr<faMesh> regionMeshPtr_;
-
         //- Model coefficients dictionary
         dictionary coeffs_;
 
+        //- Pointer to the region mesh database
+        autoPtr<faMesh> regionMeshPtr_;
+
         //- Dictionary of output properties
         autoPtr<IOdictionary> outputPropertiesPtr_;
 
@@ -224,6 +238,15 @@ public:
         //- Return mapping between surface and volume fields
         const volSurfaceMapping& vsm() const;
 
+        //- The suffix hint for automatic model variable names
+        const word& suffixHint() const noexcept { return suffixHint_; }
+
+        //- Return the concatenation of \p base and the suffix hint
+        inline word suffixed(const char* base) const;
+
+        //- Return the concatenation of \p base and the suffix hint
+        inline word suffixed(const std::string& base) const;
+
 
     // Evolution
 

src/regionFaModels/regionFaModel/regionFaModelI.H

@@ -27,7 +27,7 @@ License
 
 #include "regionFaModel.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 inline const Foam::faMesh& Foam::regionModels::regionFaModel::regionMesh() const
 {
@@ -106,4 +106,18 @@ inline bool Foam::regionModels::regionFaModel::isRegionPatch
 }
 
 
+inline Foam::word
+Foam::regionModels::regionFaModel::suffixed(const char* base) const
+{
+    return word(base + suffixHint_);
+}
+
+
+inline Foam::word
+Foam::regionModels::regionFaModel::suffixed(const std::string& base) const
+{
+    return word(base + suffixHint_);
+}
+
+
 // ************************************************************************* //

src/regionFaModels/thermalShell/thermalShell.C

@@ -123,14 +123,29 @@ thermalShell::thermalShell
 )
 :
     thermalShellModel(modelType, mesh, dict),
+    hName_(dict.getOrDefault<word>("h", suffixed("hs"))),
+    qsName_(dict.getOrDefault<word>("qs", suffixed("qs"))),
+    qrName_(dict.getOrDefault<word>("qr", "none")),
     nNonOrthCorr_(1),
     // Only need/want thermal solid properties
     thermo_(dict.subDict("thermo"), solidProperties::THERMAL),
+    h_
+    (
+        IOobject
+        (
+            hName_,
+            regionMesh().time().timeName(),
+            regionMesh().thisDb(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        regionMesh()
+    ),
     qs_
     (
         IOobject
         (
-            "qs_" + regionName_,
+            qsName_,
             regionMesh().time().timeName(),
             regionMesh().thisDb(),
             IOobject::READ_IF_PRESENT,
@@ -139,19 +154,6 @@ thermalShell::thermalShell
         regionMesh(),
         dimensionedScalar(dimPower/dimArea, Zero)
     ),
-    h_
-    (
-        IOobject
-        (
-            "h_" + regionName_,
-            regionMesh().time().timeName(),
-            regionMesh().thisDb(),
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        regionMesh()
-    ),
-    qrName_(dict.getOrDefault<word>("qr", "none")),
     thickness_(dict.getOrDefault<scalar>("thickness", 0))
 {
     init(dict);

src/regionFaModels/thermalShell/thermalShell.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -56,11 +56,28 @@ Usage
 
     where the entries mean:
     \table
-      Property   | Description                      | Type  | Reqd   | Deflt
-      thermalShellModel | Type name: thermalShell   | word  | yes    | -
-      thermo     | Solid thermal properties         | dictionary | yes | -
-      qr         | Name of radiative heat flux field | word | no     | none
-      thickness  | Uniform film thickness [m]       | scalar | choice | -
+      Property   | Description                    | Type  | Reqd   | Deflt
+      thermalShellModel | Type name: thermalShell | word  | yes    | -
+      h          | Name of thickness field        | word  | no     | hs (suffix)
+      qs         | Name of source field           | word  | no     | qs (suffix)
+      qr         | Name of radiative heat flux field | word | no   | none
+      thermo     | Solid thermal properties       | dict | yes | -
+      thickness  | Uniform shell thickness [m]    | scalar | choice | -
+    \endtable
+
+    Fields/variables used:
+    \table
+      Property  | Description               | Type  | Deflt
+      h         | Thickness                 | shell | hs (suffix)
+      qs        | Source field              | shell | qs (suffix)
+      qr        | Radiative heat flux field | volume | none
+    \endtable
+
+    Note the following naming changes from 2056 and earlier:
+    \table
+      Keyword  | Description               | Keyword (old) | Deflt (old)
+      h        | Temperature (volume)      | h    | "h_" + regionName
+      qs       | Temperature (shell)       | qs   | "qs_" + regionName
     \endtable
 
     The inherited entries are elaborated in:
@@ -68,7 +85,6 @@ Usage
 
 SourceFiles
     thermalShell.C
-    thermalShellI.H
 
 \*---------------------------------------------------------------------------*/
 
@@ -108,31 +124,38 @@ protected:
 
     // Protected Data
 
-        // Solution parameters
+        //- Name of shell thickness [height] field (default: "hs" + suffix)
+        const word hName_;
 
-            //- Number of non orthogonal correctors
-            label nNonOrthCorr_;
+        //- Name of surface energy source (default: "qs" + suffix)
+        const word qsName_;
+
+        //- Name of (volume) radiative flux field (default: none)
+        const word qrName_;
 
 
-        // Thermo properties
+    // Solution Parameters
 
-            //- Solid properties
-            solidProperties thermo_;
+        //- Number of non orthogonal correctors
+        label nNonOrthCorr_;
 
 
-        // Source term fields
+    // Thermo properties
 
-            //- External surface energy source [J/m2/s]
-            areaScalarField qs_;
+        //- Solid properties
+        solidProperties thermo_;
 
-            //- Film thickness [m]
-            areaScalarField h_;
 
-            //- Name of the primary region radiative flux
-            const word qrName_;
+    // Source term fields
 
-            //- Uniform film thickness [m]
-            scalar thickness_;
+        //- Shell thickness field [m]
+        areaScalarField h_;
+
+        //- External surface energy source [J/m2/s]
+        areaScalarField qs_;
+
+        //- Uniform shell thickness [m]
+        scalar thickness_;
 
 
     // Protected Member Functions
@@ -174,7 +197,7 @@ public:
 
         // Fields
 
-            //- Return the film specific heat capacity [J/kg/K]
+            //- Return the shell specific heat capacity [J/kg/K]
             const tmp<areaScalarField> Cp() const;
 
             //- Return density [kg/m3]

src/regionFaModels/thermalShellModel/thermalShellModel.C

@@ -50,13 +50,14 @@ thermalShellModel::thermalShellModel
 )
 :
     regionFaModel(mesh, "thermalShell", modelType, dict, true),
-    TName_(dict.get<word>("T")),
-    Tp_(mesh.lookupObject<volScalarField>(TName_)),
+    TName_(dict.getOrDefault<word>("T", suffixed("Ts"))),
+    TprimaryName_(dict.getOrDefault<word>("Tprimary", "T")),
+    Tp_(mesh.lookupObject<volScalarField>(TprimaryName_)),
     T_
     (
         IOobject
         (
-            "Ts_" + regionName_,
+            TName_,
             regionMesh().time().timeName(),
             regionMesh().thisDb(),
             IOobject::MUST_READ,
@@ -71,8 +72,8 @@ thermalShellModel::thermalShellModel
 {
     if (faOptions_.optionList::empty())
     {
-        Info<< "No finite area options present for area : "
-            << polyMesh::regionName(regionFaModel::areaName()) << endl;
+        Info<< "No finite-area options present for area:"
+            << regionFaModel::areaName() << endl;
     }
 }
 

src/regionFaModels/thermalShellModel/thermalShellModel.H

@@ -24,7 +24,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::regionModels::thermalShellModels::thermalShellModel
+    Foam::regionModels::thermalShellModel
 
 Description
     Intermediate class for thermal-shell finite-area models.
@@ -34,12 +34,12 @@ Usage
     \verbatim
     <patchName>
     {
-        // Mandatory entries
-        T                     <word>;
-
         // Optional entries
         thermalShellModel     <word>;
 
+        T                     <word>;
+        Tprimary              <word>;
+
         // Inherited entries
         ...
     }
@@ -47,12 +47,27 @@ Usage
 
     where the entries mean:
     \table
-      Property   | Description                       | Type  | Reqd | Deflt
-      T          | Name of operand temperature field | word  | yes  | -
-      thermalShellModel | Name of thermal-shell model <!--
+      Property  | Description                       | Type  | Reqd | Deflt
+      T         | Name of operand temperature field | word  | no   | Ts (suffix)
+      Tprimary  | Name of primary temperature field | word  | no   | T
+      thermalShellModel | Name of thermalShellModel thermal-shell model <!--
                         -->                      | word | choice | -
     \endtable
 
+    Fields/variables used:
+    \table
+      Property  | Description               | Type   | Deflt
+      T         | Temperature               | shell  | Ts (suffix)
+      Tprimary  | Temperature               | volume | T
+    \endtable
+
+    \b BREAKING Naming changes from 2056 and earlier:
+    \table
+      Keyword  | Description               | Keyword (old) | Deflt (old)
+      T        | Temperature (shell)       | -   | "Ts_" + regionName
+      Tprimary | Temperature (volume)      | T   | -
+    \endtable
+
     The inherited entries are elaborated in:
       - \link regionFaModel.H \endlink
 
@@ -89,8 +104,11 @@ protected:
 
     // Protected Data
 
-        //- Name of the temperature field
-        word TName_;
+        //- Name of shell temperature field (default: "Ts" + suffix)
+        const word TName_;
+
+        //- Name of volume temperature field (default: "T")
+        const word TprimaryName_;
 
         //- Primary (volume) region temperature
         const volScalarField& Tp_;

src/regionFaModels/vibrationShellModel/vibrationShellModel.C

@@ -53,7 +53,7 @@ vibrationShellModel::vibrationShellModel
     (
         IOobject
         (
-            "ws_" + regionName_,
+            dict.getOrDefault<word>("ws", suffixed("ws")),
             regionMesh().time().timeName(),
             regionMesh().thisDb(),
             IOobject::MUST_READ,
@@ -65,7 +65,7 @@ vibrationShellModel::vibrationShellModel
     (
         IOobject
         (
-            "as_" + regionName_,
+            dict.getOrDefault<word>("as", suffixed("as")),
             regionMesh().time().timeName(),
             regionMesh().thisDb(),
             IOobject::NO_READ,

src/regionFaModels/vibrationShellModel/vibrationShellModel.H

@@ -24,7 +24,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::regionModels::thermalShellModels::vibrationShellModel
+    Foam::regionModels::vibrationShellModel
 
 Description
     Intermediate class for vibration-shell finite-area models.
@@ -52,9 +52,19 @@ Usage
     \table
       Property   | Description                            | Type  | Reqd | Deflt
       vibrationShellModel | Name of vibration-shell model | word | yes | -
-      p          | Name of the coupled field in the primary <!--
-                 --> region                               | word | yes | -
-      solid      | Solid properties                 | dictionary | yes | -
+      p          | Name of primary pressure field         | word | yes | -
+      solid      | Solid properties                       | dict | yes | -
+      as         | Name of shell acceleration field    | word | no | as (suffix)
+      ps         | Name of pressure (source) on shell  | word | no | ps (suffix)
+      ws         | Name of shell displacement field    | word | no | ws (suffix)
+    \endtable
+
+    Naming changes from 2056 and earlier:
+    \table
+      Keyword  | Description               | Keyword (old) | Deflt (old)
+      as       | Acceleration (shell)      | - | "as_" + regionName
+      ps       | Pressure on shell         | - | "ps_" + regionName
+      ws       | Displacement (shell)      | - | "ws_" + regionName
     \endtable
 
     The inherited entries are elaborated in:
@@ -93,19 +103,19 @@ protected:
 
     // Protected Data
 
-        //- Shell displacement
+        //- Shell displacement [m]
         areaScalarField w_;
 
-        //- Shell acceleration
+        //- Shell acceleration [m/s2]
         areaScalarField a_;
 
         //- Solid properties
         solidProperties solid_;
 
-        //- Name of the coupled field in the primary region
+        //- Name of primary region acoustic pressure field
         word pName_;
 
-        //- Primary region acoustic pressure
+        //- Primary region acoustic pressure [Pa]
         const volScalarField& pa_;
 
         //- Reference to faOptions
@@ -167,26 +177,26 @@ public:
 
     // Member Functions
 
-       //- Return primary region pa
-       const volScalarField& pa() const noexcept { return pa_; }
+        //- Return the primary region presssure
+        const volScalarField& pa() const noexcept { return pa_; }
 
-       //- Return shell displacement
-       const areaScalarField& w() const noexcept { return w_; }
+        //- Return the shell displacement
+        const areaScalarField& w() const noexcept { return w_; }
 
-       //- Return shell acceleration
-       const areaScalarField& a() const noexcept { return a_; }
+        //- Return the shell acceleration
+        const areaScalarField& a() const noexcept { return a_; }
 
-       //- Return faOptions
-       Foam::fa::options& faOptions() noexcept { return faOptions_; }
+        //- Return faOptions
+        Foam::fa::options& faOptions() noexcept { return faOptions_; }
 
-       //- Return solid properties
-       const solidProperties& solid() const noexcept { return solid_; }
+        //- Return solid properties
+        const solidProperties& solid() const noexcept { return solid_; }
 
 
    // Evolution
 
-       //- Pre-evolve region
-       virtual void preEvolveRegion();
+        //- Pre-evolve region
+        virtual void preEvolveRegion();
 };
 
 

tutorials/compressible/acousticFoam/obliqueAirJet/main/0.orig/finite-area/hs

@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       areaScalarField;
-    object      h_ceilingShell;
+    object      hs;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

tutorials/compressible/acousticFoam/obliqueAirJet/main/0.orig/finite-area/ws

@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       areaScalarField;
-    object      wS;
+    object      ws;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

tutorials/compressible/acousticFoam/obliqueAirJet/main/0.orig/pa

@@ -41,7 +41,7 @@ boundaryField
             nu      0.22;
         }
 
-        region                vibrationShell;
+        region                vibration;
         vibrationShellModel   KirchhoffShell;
 
         f0             0.04;

tutorials/compressible/acousticFoam/obliqueAirJet/main/system/finite-area/faOptions

@@ -17,12 +17,13 @@ FoamFile
 pressure
 {
     type            externalFileSource;
-    fieldName       ws_vibrationShell;
+    //fieldName       ws_shell;
+    fieldName       ws;
     tableName       p;
     active          true;
     timeStart       0.001;
     duration        0.03;
-    region          vibrationShell;
+    region          vibration;
     selectionMode   all;
 
     sampleFormat    ensight;