WIP: ENH: evapotranspirationHeatTransfer: New atmospheric-model fvOption

WIP: The units of the main governing eq are inconsistent - see EP1950

Applies sources on temperature ('T' - incompressible)
or energy ('h'/'e' - compressible) equation to incorporate
evapotranspiration heat-transfer effects from the specified
plant canopy. Heat transfer is usually calculated based on
empirical relations between plants and solar radiation.

Two submodels to incorporate heat transfer effects
- 'tree': specified tree canopy - uses empirical relations between solar
radiation and evapotranspiration.
- 'grass': specified grass canopy - uses Pemnan-Monteith Equation model.

applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C

@@ -145,6 +145,7 @@ Usage
 #include "hexRef8Data.H"
 #include "regionProperties.H"
 #include "polyMeshTools.H"
+#include "subsetAdjacency.H"
 
 using namespace Foam;
 
@@ -194,21 +195,56 @@ tmp<volScalarField> createScalarField
 }
 
 
-// Calculate band of matrix
-label getBand(const labelList& owner, const labelList& neighbour)
+// Calculate band of mesh
+// label getBand(const labelUList& owner, const labelUList& neighbour)
+// {
+//     label bandwidth = 0;
+//
+//     forAll(neighbour, facei)
+//     {
+//         const label width = neighbour[facei] - owner[facei];
+//
+//         if (bandwidth < width)
+//         {
+//             bandwidth = width;
+//         }
+//     }
+//     return bandwidth;
+// }
+
+
+// Calculate band and profile of matrix. Profile is scalar to avoid overflow
+Tuple2<label, scalar> getBand
+(
+    const CompactListList<label>& mat
+)
 {
-    label band = 0;
+    Tuple2<label, scalar> metrics(0, 0);
 
-    forAll(neighbour, facei)
+    auto& bandwidth = metrics.first();
+    auto& profile = metrics.second();
+
+    forAll(mat, celli)
     {
-        label diff = neighbour[facei] - owner[facei];
+        const auto& neighbours = mat[celli];
 
-        if (diff > band)
+        const label nNbr = neighbours.size();
+
+        if (nNbr)
         {
-            band = diff;
+            // Max distance
+            const label width = (neighbours[nNbr-1] - celli);
+
+            if (bandwidth < width)
+            {
+                bandwidth = width;
+            }
+
+            profile += scalar(width);
         }
     }
-    return band;
+
+    return metrics;
 }
 
 
@@ -217,27 +253,35 @@ void getBand
 (
     const bool calculateIntersect,
     const label nCells,
-    const labelList& owner,
-    const labelList& neighbour,
+    const labelUList& owner,
+    const labelUList& neighbour,
     label& bandwidth,
     scalar& profile,            // scalar to avoid overflow
     scalar& sumSqrIntersect     // scalar to avoid overflow
 )
 {
     labelList cellBandwidth(nCells, Foam::zero{});
-    scalarField nIntersect(nCells, Foam::zero{});
+
+    bandwidth = 0;
 
     forAll(neighbour, facei)
     {
-        label own = owner[facei];
-        label nei = neighbour[facei];
+        const label own = owner[facei];
+        const label nei = neighbour[facei];
 
         // Note: mag not necessary for correct (upper-triangular) ordering.
-        label diff = nei-own;
-        cellBandwidth[nei] = max(cellBandwidth[nei], diff);
-    }
+        const label width = nei - own;
 
-    bandwidth = max(cellBandwidth);
+        if (cellBandwidth[nei] < width)
+        {
+            cellBandwidth[nei] = width;
+
+            if (bandwidth < width)
+            {
+                bandwidth = width;
+            }
+        }
+    }
 
     // Do not use field algebra because of conversion label to scalar
     profile = 0;
@@ -246,14 +290,16 @@ void getBand
         profile += scalar(width);
     }
 
-    sumSqrIntersect = 0.0;
+    sumSqrIntersect = 0;
     if (calculateIntersect)
     {
+        scalarField nIntersect(nCells, Foam::zero{});
+
         forAll(nIntersect, celli)
         {
             for (label colI = celli-cellBandwidth[celli]; colI <= celli; colI++)
             {
-                nIntersect[colI] += 1.0;
+                nIntersect[colI] += scalar(1);
             }
         }
 
@@ -675,9 +721,8 @@ CompactListList<label> regionRenumber
 
         forAll(regionCellOrder, regioni)
         {
-            // Info<< "    region " << regioni
-            //     << " starts at " << regionCellOrder.localStart(regioni)
-            //     << nl;
+            // Info<< "    region " << regioni << " starts at "
+            //     << regionCellOrder.localStart(regioni) << nl;
 
             // No parallel communication
             const bool oldParRun = UPstream::parRun(false);
@@ -703,27 +748,42 @@ CompactListList<label> regionRenumber
     {
         timer.resetTimeIncrement();
 
+        // Create adjacency matrix of the full mesh and subset subsequently.
+        // This is more efficient than creating adjacency matrices of
+        // sub-meshes.
+
+        // No parallel communication
+        const bool oldParRun = UPstream::parRun(false);
+
+        // The local connectivity of the full (non-subsetted) mesh
+        CompactListList<label> meshCellCells;
+        globalMeshData::calcCellCells(mesh, meshCellCells);
+        UPstream::parRun(oldParRun);  // Restore parallel state
+
+        timings[TimingType::CELL_CELLS] += timer.timeIncrement();
+
+        // For the respective subMesh selections
+        bitSet subsetCells(mesh.nCells());
+
         forAll(regionCellOrder, regioni)
         {
-            // Info<< "    region " << regioni
-            //     << " starts at " << regionCellOrder.localStart(regioni)
-            //     << nl;
+            // Info<< "    region " << regioni << " starts at "
+            //     << regionCellOrder.localStart(regioni) << nl;
+
+            subsetCells = false;
+            subsetCells.set(regionCellOrder[regioni]);
+
+            // Connectivity of local sub-mesh
+            labelList cellMap;
+            CompactListList<label> subCellCells =
+                subsetAdjacency(subsetCells, meshCellCells, cellMap);
+
+            timings[TimingType::CELL_CELLS] += timer.timeIncrement();
 
             // No parallel communication
             const bool oldParRun = UPstream::parRun(false);
 
-            // Connectivity of local sub-mesh
-            CompactListList<label> cellCells;
-            labelList cellMap = globalMeshData::calcCellCells
-            (
-                mesh,
-                regionCellOrder[regioni],
-                cellCells
-            );
-
-            timings[TimingType::CELL_CELLS] += timer.timeIncrement();
-
-            labelList subCellOrder = method.renumber(cellCells);
+            labelList subCellOrder = method.renumber(subCellCells);
 
             UPstream::parRun(oldParRun);  // Restore parallel state
 
@@ -835,10 +895,10 @@ int main(int argc, char *argv[])
     const bool dryrun = args.dryRun();
 
     const bool readDict = args.found("dict");
-    const bool doFrontWidth = args.found("frontWidth");
+    const bool doDecompose = args.found("decompose");
     const bool overwrite = args.found("overwrite");
     const bool doFields = !args.found("no-fields");
-    const bool doDecompose = args.found("decompose");
+    const bool doFrontWidth = args.found("frontWidth") && !doDecompose;
 
     word renumberMethodName;
     args.readIfPresent("renumber-method", renumberMethodName);
@@ -846,8 +906,7 @@ int main(int argc, char *argv[])
     if (doDecompose && UPstream::parRun())
     {
         FatalErrorIn(args.executable())
-            << "Cannot use -decompose option in parallel"
-            << " ... giving up" << nl
+            << "Cannot use -decompose option in parallel ... giving up" << nl
             << exit(FatalError);
     }
 
@@ -908,21 +967,21 @@ int main(int argc, char *argv[])
 
         reduce(band, maxOp<label>());
         reduce(profile, sumOp<scalar>());
-        reduce(sumSqrIntersect, sumOp<scalar>());
-
-        scalar rmsFrontwidth = Foam::sqrt
-        (
-            sumSqrIntersect/mesh.globalData().nTotalCells()
-        );
 
         Info<< "Mesh " << mesh.name()
             << " size: " << mesh.globalData().nTotalCells() << nl
-            << "Before renumbering :" << nl
+            << "Before renumbering" << nl
             << "    band           : " << band << nl
             << "    profile        : " << profile << nl;
 
         if (doFrontWidth)
         {
+            reduce(sumSqrIntersect, sumOp<scalar>());
+            scalar rmsFrontwidth = Foam::sqrt
+            (
+                sumSqrIntersect/mesh.globalData().nTotalCells()
+            );
+
             Info<< "    rms frontwidth : " << rmsFrontwidth << nl;
         }
 
@@ -1091,10 +1150,7 @@ int main(int argc, char *argv[])
         );
 
 
-        // List of objects read from time directory
-        // List of stored objects to clear from mesh
-
-        IOobjectList objects;
+        // List of stored objects to clear from mesh (after reading)
         DynamicList<regIOobject*> storedObjects;
 
         if (!dryrun && doFields)
@@ -1103,55 +1159,60 @@ int main(int argc, char *argv[])
 
             timer.resetTimeIncrement();
 
-            objects = IOobjectList(mesh, runTime.timeName());
+            IOobjectList objects(mesh, runTime.timeName());
             storedObjects.reserve(objects.size());
 
             const predicates::always nameMatcher;
 
             // Read GeometricFields
 
-            #undef  ReadFields
-            #define ReadFields(FieldType)                                     \
+            #undef  doLocalCode
+            #define doLocalCode(FieldType)                                    \
             readFields<FieldType>(mesh, objects, nameMatcher, storedObjects);
 
             // Read volume fields
-            ReadFields(volScalarField);
-            ReadFields(volVectorField);
-            ReadFields(volSphericalTensorField);
-            ReadFields(volSymmTensorField);
-            ReadFields(volTensorField);
+            doLocalCode(volScalarField);
+            doLocalCode(volVectorField);
+            doLocalCode(volSphericalTensorField);
+            doLocalCode(volSymmTensorField);
+            doLocalCode(volTensorField);
 
             // Read internal fields
-            ReadFields(volScalarField::Internal);
-            ReadFields(volVectorField::Internal);
-            ReadFields(volSphericalTensorField::Internal);
-            ReadFields(volSymmTensorField::Internal);
-            ReadFields(volTensorField::Internal);
+            doLocalCode(volScalarField::Internal);
+            doLocalCode(volVectorField::Internal);
+            doLocalCode(volSphericalTensorField::Internal);
+            doLocalCode(volSymmTensorField::Internal);
+            doLocalCode(volTensorField::Internal);
 
             // Read surface fields
-            ReadFields(surfaceScalarField);
-            ReadFields(surfaceVectorField);
-            ReadFields(surfaceSphericalTensorField);
-            ReadFields(surfaceSymmTensorField);
-            ReadFields(surfaceTensorField);
+            doLocalCode(surfaceScalarField);
+            doLocalCode(surfaceVectorField);
+            doLocalCode(surfaceSphericalTensorField);
+            doLocalCode(surfaceSymmTensorField);
+            doLocalCode(surfaceTensorField);
 
 
             // Read point fields
             const pointMesh& pMesh = pointMesh::New(mesh);
-            #undef  ReadPointFields
-            #define ReadPointFields(FieldType)                                \
+            #undef  doLocalCode
+            #define doLocalCode(FieldType)                                    \
             readFields<FieldType>(pMesh, objects, nameMatcher, storedObjects);
 
-            ReadPointFields(pointScalarField);
-            ReadPointFields(pointVectorField);
-            ReadPointFields(pointSphericalTensorField);
-            ReadPointFields(pointSymmTensorField);
-            ReadPointFields(pointTensorField);
+            doLocalCode(pointScalarField);
+            doLocalCode(pointVectorField);
+            doLocalCode(pointSphericalTensorField);
+            doLocalCode(pointSymmTensorField);
+            doLocalCode(pointTensorField);
 
-            #undef ReadFields
-            #undef ReadPointFields
+            #undef doLocalCode
 
             timings[TimingType::READ_FIELDS] += timer.timeIncrement();
+
+            // Write loaded fields when mesh.write() is called
+            for (auto* fldptr : storedObjects)
+            {
+                fldptr->writeOpt(IOobject::AUTO_WRITE);
+            }
         }
 
 
@@ -1224,7 +1285,13 @@ int main(int argc, char *argv[])
 
 
             CompactListList<label> regionCellOrder =
-                regionRenumber(renumberPtr(), mesh, cellToRegion);
+                regionRenumber
+                (
+                    renumberPtr(),
+                    mesh,
+                    cellToRegion,
+                    decomposePtr().nDomains()
+                );
 
             cellOrder = regionCellOrder.values();
 
@@ -1586,19 +1653,25 @@ int main(int argc, char *argv[])
             );
             reduce(band, maxOp<label>());
             reduce(profile, sumOp<scalar>());
-            reduce(sumSqrIntersect, sumOp<scalar>());
 
-            scalar rmsFrontwidth = Foam::sqrt
-            (
-                sumSqrIntersect/mesh.globalData().nTotalCells()
-            );
+            Info<< "After renumbering";
+            if (doDecompose)
+            {
+                Info<< " [values are misleading with -decompose option]";
+            }
 
-            Info<< "After renumbering :" << nl
+            Info<< nl
                 << "    band           : " << band << nl
                 << "    profile        : " << profile << nl;
 
             if (doFrontWidth)
             {
+                reduce(sumSqrIntersect, sumOp<scalar>());
+                scalar rmsFrontwidth = Foam::sqrt
+                (
+                    sumSqrIntersect/mesh.globalData().nTotalCells()
+                );
+
                 Info<< "    rms frontwidth : " << rmsFrontwidth << nl;
             }
 

applications/utilities/mesh/manipulation/renumberMesh/subsetAdjacency.H

@@ -0,0 +1,168 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2024 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+   This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
+
+Description
+    Subsetting of an adjacency matrix (as CompactListList).
+    Can be relocated elsewhere.
+
+\*---------------------------------------------------------------------------*/
+
+#include "CompactListList.H"
+#include "bitSet.H"
+#include "ListOps.H"
+#include "Map.H"
+
+namespace Foam
+{
+
+// Perform a subset of the adjacency matrix
+CompactListList<label> subsetAdjacency
+(
+    const bitSet& select,   // could also be labelHashSet
+    const CompactListList<label>& input,
+    labelList& subMap
+)
+{
+    // Corresponds to cellMap etc (the original selection)
+    subMap = select.sortedToc();
+
+    // Ensure that the subMap corresponds to a valid subset
+    {
+        label validSize = 0;
+
+        const label nTotal = input.size();
+
+        forAllReverse(subMap, i)
+        {
+            if (subMap[i] < nTotal)
+            {
+                validSize = i + 1;
+                break;
+            }
+        }
+
+        subMap.resize(validSize);
+    }
+
+
+    // Assumed to be sparse - use Map for reverse lookup
+    const Map<label> reverseMap(invertToMap(subMap));
+
+
+    // Pass 1: determine the selected sub-sizes
+    labelList sizes(subMap.size(), Foam::zero{});
+
+    forAll(subMap, idx)
+    {
+        for (const label nbr : input[subMap[idx]])
+        {
+            if
+            (
+                select.test(nbr)
+             && reverseMap.contains(nbr)  // extra consistency (paranoid)
+            )
+            {
+                ++sizes[idx];
+            }
+        }
+    }
+
+
+    CompactListList<label> output(sizes);
+
+    // Reuse sizes as output offset into output.values()
+    sizes = labelList::subList(output.offsets(), output.size());
+    labelList& values = output.values();
+
+
+    // Pass 2: extract sub-adjacent matrix
+
+    label newNbr = -1;
+
+    forAll(subMap, idx)
+    {
+        for (const label nbr : input[subMap[idx]])
+        {
+            if
+            (
+                select.test(nbr)
+             && (newNbr = reverseMap.lookup(nbr, -1)) >= 0
+            )
+            {
+                values[sizes[idx]++] = newNbr;
+            }
+        }
+    }
+
+    return output;
+}
+
+
+// Perform a subset of the adjacency matrix
+CompactListList<label> subsetAdjacency
+(
+    const labelRange& slice,
+    const CompactListList<label>& input,
+    labelList& subMap
+)
+{
+    // Ensure that the selection corresponds to a valid subset
+    const labelRange select = slice.subset0(input.size());
+
+    // Corresponds to cellMap etc (the original selection)
+    subMap = Foam::identity(select);
+
+
+    // Pass 1: determine the selected sub-sizes
+    labelList sizes(subMap.size(), Foam::zero{});
+
+    forAll(subMap, idx)
+    {
+        for (const label nbr : input[subMap[idx]])
+        {
+            if (select.contains(nbr))
+            {
+                ++sizes[idx];
+            }
+        }
+    }
+
+
+    CompactListList<label> output(sizes);
+
+    // Reuse sizes as output offset into output.values()
+    sizes = labelList::subList(output.offsets(), output.size());
+    labelList& values = output.values();
+
+
+    // Pass 2: extract sub-adjacent matrix
+
+    const label localOffset = select.start();
+
+    forAll(subMap, idx)
+    {
+        for (const label nbr : input[subMap[idx]])
+        {
+            if (select.contains(nbr))
+            {
+                values[sizes[idx]++] = nbr - localOffset;
+            }
+        }
+    }
+
+    return output;
+}
+
+} // End namespace Foam
+
+
+// ************************************************************************* //

src/OpenFOAM/containers/Bits/PackedList/PackedList.H

@@ -83,8 +83,9 @@ See also
     Foam::DynamicList
 
 SourceFiles
-    PackedListI.H
     PackedList.C
+    PackedListCore.C
+    PackedListI.H
     PackedListIO.C
 
 \*---------------------------------------------------------------------------*/
@@ -92,11 +93,11 @@ SourceFiles
 #ifndef Foam_PackedList_H
 #define Foam_PackedList_H
 
+#include "className.H"
 #include "BitOps.H"
 #include "labelList.H"
 #include "IndirectListBase.H"
 #include "InfoProxy.H"
-#include "PackedListCore.H"
 
 #include <type_traits>
 
@@ -119,6 +120,23 @@ template<unsigned Width>
 Ostream& operator<<(Ostream& os, const InfoProxy<PackedList<Width>>& info);
 
 
+/*---------------------------------------------------------------------------*\
+                   Class Detail::PackedListCore Declaration
+\*---------------------------------------------------------------------------*/
+
+namespace Detail
+{
+
+//- Template-invariant parts for PackedList
+struct PackedListCore
+{
+    //- Define template name
+    ClassNameNoDebug("PackedList");
+};
+
+} // End namespace Detail
+
+
 /*---------------------------------------------------------------------------*\
                          Class PackedList Declaration
 \*---------------------------------------------------------------------------*/

src/OpenFOAM/containers/Bits/PackedList/PackedListCore.C

@@ -25,7 +25,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "PackedListCore.H"
+#include "PackedList.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 

src/OpenFOAM/containers/Bits/PackedList/PackedListCore.H

@@ -1,69 +1 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2018 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::Detail::PackedListCore
-
-Description
-    Implementation of template-invariant details for Foam::PackedList
-
-SourceFiles
-    PackedListCore.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef PackedListCore_H
-#define PackedListCore_H
-
-#include "className.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace Detail
-{
-
-/*---------------------------------------------------------------------------*\
-                   Class Detail::PackedListCore Declaration
-\*---------------------------------------------------------------------------*/
-
-//- Template-invariant parts for PackedList
-struct PackedListCore
-{
-    //- Define template name
-    ClassNameNoDebug("PackedList");
-};
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Detail
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
+#warning File removed - left for old dependency check only

src/OpenFOAM/global/fileOperations/uncollatedFileOperation/uncollatedFileOperation.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2017 OpenFOAM Foundation
-    Copyright (C) 2020-2023 OpenCFD Ltd.
+    Copyright (C) 2020-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -716,6 +716,14 @@ bool Foam::fileOperations::uncollatedFileOperation::read
         const bool oldMasterOnly = regIOobject::masterOnlyReading;
         regIOobject::masterOnlyReading = masterOnly;
 
+        const bool oldParRun = UPstream::parRun();
+        if (masterOnly)
+        {
+            // Reading on master only.
+            // Avoid side effects from io.readStream below.
+            UPstream::parRun(false);
+        }
+
         // Read file
         ok = io.readData(io.readStream(typeName));
         io.close();
@@ -723,6 +731,7 @@ bool Foam::fileOperations::uncollatedFileOperation::read
         // Restore flags
         io.globalObject(oldGlobal);
         regIOobject::masterOnlyReading = oldMasterOnly;
+        UPstream::parRun(oldParRun);
 
         if (debug)
         {

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshDataTopology.C

@@ -84,7 +84,10 @@ static void calcCellCellsImpl
             globalNeighbour[bfacei] = val;
         }
     }
-    syncTools::swapBoundaryFaceList(mesh, globalNeighbour);
+
+    // Swap boundary neighbour information:
+    // - cyclics and (optionally) processor
+    syncTools::swapBoundaryFaceList(mesh, globalNeighbour, parallel);
 
 
     // Count number of faces (internal + coupled)

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeBaseTemplates.C

@@ -349,7 +349,10 @@ void Foam::mapDistributeBase::receive
             {
                 FatalErrorInFunction
                     << "From processor " << proci
-                    << " : unallocated receive field" << nl
+                    << " : unallocated receive field."
+                    << " Expected size " << map.size()
+                    << " on comm " << comm
+                    << " with procs " << UPstream::nProcs(comm) << nl
                     << exit(FatalError);
             }
         }

src/OpenFOAM/meshes/polyMesh/syncTools/syncTools.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -34,34 +34,34 @@ void Foam::syncTools::swapBoundaryCellPositions
 (
     const polyMesh& mesh,
     const UList<point>& cellData,
-    List<point>& neighbourCellData
+    List<point>& neighbourCellData,
+    const bool parRun
 )
 {
     if (cellData.size() != mesh.nCells())
     {
         FatalErrorInFunction
-            << "Number of cell values " << cellData.size()
-            << " is not equal to the number of cells in the mesh "
-            << mesh.nCells() << abort(FatalError);
+            << "Number of values " << cellData.size()
+            << " != number of cells " << mesh.nCells() << nl
+            << abort(FatalError);
     }
 
-    const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
     neighbourCellData.resize(mesh.nBoundaryFaces());
 
-    for (const polyPatch& pp : patches)
+    for (const polyPatch& pp : mesh.boundaryMesh())
     {
-        label bFacei = pp.start()-mesh.nInternalFaces();
+        const auto& faceCells = pp.faceCells();
 
-        const labelUList& faceCells = pp.faceCells();
-
-        for (const label celli : faceCells)
-        {
-            neighbourCellData[bFacei] = cellData[celli];
-            ++bFacei;
-        }
+        // ie, boundarySlice() = patchInternalList()
+        SubList<point>
+        (
+            neighbourCellData,
+            faceCells.size(),
+            pp.offset()
+        ) = UIndirectList<point>(cellData, faceCells);
     }
-    syncTools::swapBoundaryFacePositions(mesh, neighbourCellData);
+
+    syncTools::swapBoundaryFacePositions(mesh, neighbourCellData, parRun);
 }
 
 
@@ -127,9 +127,7 @@ Foam::bitSet Foam::syncTools::getMasterFaces(const polyMesh& mesh)
 {
     bitSet isMaster(mesh.nFaces(), true);
 
-    const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-    for (const polyPatch& pp : patches)
+    for (const polyPatch& pp : mesh.boundaryMesh())
     {
         if (pp.coupled())
         {
@@ -151,9 +149,7 @@ Foam::bitSet Foam::syncTools::getInternalOrMasterFaces
 {
     bitSet isMaster(mesh.nFaces(), true);
 
-    const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-    for (const polyPatch& pp : patches)
+    for (const polyPatch& pp : mesh.boundaryMesh())
     {
         if (pp.coupled())
         {
@@ -179,9 +175,7 @@ Foam::bitSet Foam::syncTools::getInternalOrCoupledFaces
 {
     bitSet isMaster(mesh.nFaces(), true);
 
-    const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-    for (const polyPatch& pp : patches)
+    for (const polyPatch& pp : mesh.boundaryMesh())
     {
         if (!pp.coupled())
         {

src/OpenFOAM/meshes/polyMesh/syncTools/syncTools.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2021 OpenCFD Ltd.
+    Copyright (C) 2015-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -69,6 +69,7 @@ class syncTools
     // Private Member Functions
 
         //- Combine val with existing value in pointValues map at given index
+        //  No communication
         template<class T, class CombineOp>
         static void combine
         (
@@ -79,6 +80,7 @@ class syncTools
         );
 
         //- Combine val with existing value in edgeValues at edge index
+        //  No communication
         template<class T, class CombineOp>
         static void combine
         (
@@ -95,6 +97,7 @@ public:
         // Preferably use specialisations below.
 
             //- Synchronize values on selected points.
+            //  Communication if UPstream::parRun() == true.
             template<class T, class CombineOp, class TransformOp>
             static void syncPointMap
             (
@@ -105,6 +108,7 @@ public:
             );
 
             //- Synchronize values on selected edges.
+            //  Communication if UPstream::parRun() == true.
             template<class T, class CombineOp, class TransformOp>
             static void syncEdgeMap
             (
@@ -115,6 +119,7 @@ public:
             );
 
             //- Synchronize values on all mesh points.
+            //  Communication if UPstream::parRun() == true.
             template<class T, class CombineOp, class TransformOp>
             static void syncPointList
             (
@@ -126,6 +131,7 @@ public:
             );
 
             //- Synchronize values on selected mesh points.
+            //  Communication if UPstream::parRun() == true.
             template<class T, class CombineOp, class TransformOp>
             static void syncPointList
             (
@@ -138,6 +144,7 @@ public:
             );
 
             //- Synchronize values on all mesh edges.
+            //  Communication if UPstream::parRun() == true.
             template<class T, class CombineOp, class TransformOp, class FlipOp>
             static void syncEdgeList
             (
@@ -150,11 +157,12 @@ public:
             );
 
             //- Synchronize values on selected mesh edges.
+            //  Communication if UPstream::parRun() == true.
             template<class T, class CombineOp, class TransformOp, class FlipOp>
             static void syncEdgeList
             (
                 const polyMesh& mesh,
-                const labelList& meshEdges,
+                const labelUList& meshEdges,
                 List<T>& edgeValues,
                 const CombineOp& cop,
                 const T& nullValue,
@@ -170,6 +178,7 @@ public:
                 UList<T>& faceValues,
                 const CombineOp& cop,
                 const TransformOp& top,
+                //! Allow parallel communication
                 const bool parRun = UPstream::parRun()
             );
 
@@ -177,6 +186,7 @@ public:
         // Synchronise point-wise data
 
             //- Synchronize values on all mesh points.
+            //  Communication if UPstream::parRun() == true.
             template<class T, class CombineOp>
             static void syncPointList
             (
@@ -197,6 +207,7 @@ public:
             }
 
             //- Synchronize locations on all mesh points.
+            //  Communication if UPstream::parRun() == true.
             template<class CombineOp>
             static void syncPointPositions
             (
@@ -217,11 +228,12 @@ public:
             }
 
             //- Synchronize values on selected mesh points.
+            //  Communication if UPstream::parRun() == true.
             template<class T, class CombineOp>
             static void syncPointList
             (
                 const polyMesh& mesh,
-                const labelList& meshPoints,
+                const labelUList& meshPoints,
                 List<T>& pointValues,
                 const CombineOp& cop,
                 const T& nullValue
@@ -239,11 +251,12 @@ public:
             }
 
             //- Synchronize locations on selected mesh points.
+            //  Communication if UPstream::parRun() == true.
             template<class CombineOp>
             static void syncPointPositions
             (
                 const polyMesh& mesh,
-                const labelList& meshPoints,
+                const labelUList& meshPoints,
                 List<point>& positions,
                 const CombineOp& cop,
                 const point& nullValue
@@ -264,6 +277,7 @@ public:
         // Synchronise edge-wise data
 
             //- Synchronize values on all mesh edges.
+            //  Communication if UPstream::parRun() == true.
             template<class T, class CombineOp>
             static void syncEdgeList
             (
@@ -285,6 +299,7 @@ public:
             }
 
             //- Synchronize locations on all mesh edges.
+            //  Communication if UPstream::parRun() == true.
             template<class CombineOp>
             static void syncEdgePositions
             (
@@ -306,11 +321,12 @@ public:
             }
 
             //- Synchronize values on selected mesh edges.
+            //  Communication if UPstream::parRun() == true.
             template<class T, class CombineOp>
             static void syncEdgeList
             (
                 const polyMesh& mesh,
-                const labelList& meshEdges,
+                const labelUList& meshEdges,
                 List<T>& edgeValues,
                 const CombineOp& cop,
                 const T& nullValue
@@ -329,11 +345,12 @@ public:
             }
 
             //- Synchronize locations on selected mesh edges.
+            //  Communication if UPstream::parRun() == true.
             template<class CombineOp>
             static void syncEdgePositions
             (
                 const polyMesh& mesh,
-                const labelList& meshEdges,
+                const labelUList& meshEdges,
                 List<point>& positions,
                 const CombineOp& cop,
                 const point& nullValue
@@ -397,7 +414,9 @@ public:
             (
                 const polyMesh& mesh,
                 UList<T>& faceValues,
-                const CombineOp& cop
+                const CombineOp& cop,
+                //! Allow parallel communication
+                const bool parRun = UPstream::parRun()
             )
             {
                 SubList<T> bndValues
@@ -412,7 +431,8 @@ public:
                     mesh,
                     bndValues,
                     cop,
-                    mapDistribute::transform()
+                    mapDistribute::transform(),
+                    parRun
                 );
             }
 
@@ -422,7 +442,9 @@ public:
             (
                 const polyMesh& mesh,
                 UList<point>& positions,
-                const CombineOp& cop
+                const CombineOp& cop,
+                //! Allow parallel communication
+                const bool parRun = UPstream::parRun()
             )
             {
                 SubList<point> bndValues
@@ -436,7 +458,8 @@ public:
                     mesh,
                     bndValues,
                     cop,
-                    mapDistribute::transformPosition()
+                    mapDistribute::transformPosition(),
+                    parRun
                 );
             }
 
@@ -445,7 +468,9 @@ public:
             static void swapBoundaryFaceList
             (
                 const polyMesh& mesh,
-                UList<T>& faceValues
+                UList<T>& faceValues,
+                //! Allow parallel communication
+                const bool parRun = UPstream::parRun()
             )
             {
                 syncBoundaryFaceList
@@ -453,7 +478,8 @@ public:
                     mesh,
                     faceValues,
                     eqOp<T>(),
-                    mapDistribute::transform()
+                    mapDistribute::transform(),
+                    parRun
                 );
             }
 
@@ -461,7 +487,9 @@ public:
             static void swapBoundaryFacePositions
             (
                 const polyMesh& mesh,
-                UList<point>& positions
+                UList<point>& positions,
+                //! Allow parallel communication
+                const bool parRun = UPstream::parRun()
             )
             {
                 syncBoundaryFaceList
@@ -469,7 +497,8 @@ public:
                     mesh,
                     positions,
                     eqOp<point>(),
-                    mapDistribute::transformPosition()
+                    mapDistribute::transformPosition(),
+                    parRun
                 );
             }
 
@@ -478,7 +507,9 @@ public:
             static void swapFaceList
             (
                 const polyMesh& mesh,
-                UList<T>& faceValues
+                UList<T>& faceValues,
+                //! Allow parallel communication
+                const bool parRun = UPstream::parRun()
             )
             {
                 SubList<T> bndValues
@@ -492,31 +523,71 @@ public:
                     mesh,
                     bndValues,
                     eqOp<T>(),
-                    mapDistribute::transform()
+                    mapDistribute::transform(),
+                    parRun
                 );
             }
 
-            //- Swap to obtain neighbour cell values for all boundary faces
+            //- Extract and swap to obtain neighbour cell values
+            //- for all boundary faces
             template<class T>
             static void swapBoundaryCellList
             (
                 const polyMesh& mesh,
                 const UList<T>& cellData,
-                List<T>& neighbourCellData
+                List<T>& neighbourCellData,
+                //! Allow parallel communication
+                const bool parRun = UPstream::parRun()
             );
 
-            //- Swap to obtain neighbour cell positions for all boundary faces
+            //- Extract and swap to obtain neighbour cell positions
+            //- for all boundary faces
             static void swapBoundaryCellPositions
             (
                 const polyMesh& mesh,
                 const UList<point>& cellData,
-                List<point>& neighbourCellData
+                List<point>& neighbourCellData,
+                //! Allow parallel communication
+                const bool parRun = UPstream::parRun()
             );
 
 
+            //- Return neighbour cell values for all boundary faces
+            //- by swapping via boundary faces
+            template<class T>
+            FOAM_NODISCARD static List<T> swapBoundaryCellList
+            (
+                const polyMesh& mesh,
+                const UList<T>& cellData,
+                //! Allow parallel communication
+                const bool parRun = UPstream::parRun()
+            )
+            {
+                List<T> nbrCellData;
+                swapBoundaryCellList(mesh, cellData, nbrCellData, parRun);
+                return nbrCellData;
+            }
+
+            //- Return neighbour cell positions for all boundary faces
+            //- by swapping via boundary faces
+            FOAM_NODISCARD static List<point> swapBoundaryCellPositions
+            (
+                const polyMesh& mesh,
+                const UList<point>& cellData,
+                //! Allow parallel communication
+                const bool parRun = UPstream::parRun()
+            )
+            {
+                List<point> nbrCellData;
+                swapBoundaryCellPositions(mesh, cellData, nbrCellData, parRun);
+                return nbrCellData;
+            }
+
+
         // Sparse versions
 
             //- Synchronize values on selected points.
+            //  Communication if UPstream::parRun() == true.
             template<class T, class CombineOp>
             static void syncPointMap
             (
@@ -535,6 +606,7 @@ public:
             }
 
             //- Synchronize locations on selected points.
+            //  Communication if UPstream::parRun() == true.
             template<class CombineOp>
             static void syncPointPositions
             (
@@ -552,9 +624,10 @@ public:
                 );
             }
 
-            //- Synchronize values on selected edges. Edges are represented
-            //  by the two vertices that make it up so global edges never get
-            //  constructed.
+            //- Synchronize values on selected edges.
+            //- Edges are represented by the two vertices that make it up
+            //- so global edges never get constructed.
+            //  Communication if UPstream::parRun() == true.
             template<class T, class CombineOp>
             static void syncEdgeMap
             (
@@ -573,6 +646,7 @@ public:
             }
 
             //- Synchronize locations on selected edges.
+            //  Communication if UPstream::parRun() == true.
             template<class CombineOp>
             static void syncEdgePositions
             (
@@ -601,7 +675,7 @@ public:
             //      offset when accessing values.
             //  \param faceValues The face values to synchronize
             //  \param cop The combine operation
-            //  \param parRun True if this is a parallel simulation
+            //  \param parRun Allow parallel communication
             template<unsigned Width, class CombineOp>
             static void syncFaceList
             (
@@ -619,6 +693,7 @@ public:
                 const polyMesh& mesh,
                 PackedList<Width>& faceValues,
                 const CombineOp& cop,
+                //! Allow parallel communication
                 const bool parRun = UPstream::parRun()
             );
 
@@ -629,6 +704,7 @@ public:
                 const polyMesh& mesh,
                 PackedList<Width>& faceValues,
                 const CombineOp& cop,
+                //! Allow parallel communication
                 const bool parRun = UPstream::parRun()
             );
 
@@ -637,7 +713,9 @@ public:
             static void swapFaceList
             (
                 const polyMesh& mesh,
-                PackedList<Width>& faceValues
+                PackedList<Width>& faceValues,
+                //! Allow parallel communication
+                const bool parRun = UPstream::parRun()
             );
 
             //- Swap coupled boundary face values. Uses eqOp
@@ -645,7 +723,9 @@ public:
             static void swapBoundaryFaceList
             (
                 const polyMesh& mesh,
-                PackedList<Width>& faceValues
+                PackedList<Width>& faceValues,
+                //! Allow parallel communication
+                const bool parRun = UPstream::parRun()
             );
 
             template<unsigned Width, class CombineOp>

src/OpenFOAM/meshes/polyMesh/syncTools/syncToolsTemplates.C

@@ -79,6 +79,8 @@ void Foam::syncTools::combine
 }
 
 
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
 template<class T, class CombineOp, class TransformOp>
 void Foam::syncTools::syncPointMap
 (
@@ -488,7 +490,7 @@ void Foam::syncTools::syncEdgeMap
                 forAllConstIters(nbrPatchInfo, nbrIter)
                 {
                     const edge& e = nbrIter.key();
-                    const edge meshEdge(meshPts[e[0]], meshPts[e[1]]);
+                    const edge meshEdge(meshPts, e);
 
                     combine
                     (
@@ -535,7 +537,7 @@ void Foam::syncTools::syncEdgeMap
 
                 {
                     const edge& e0 = edgesA[twoEdges[0]];
-                    const edge meshEdge0(meshPtsA[e0[0]], meshPtsA[e0[1]]);
+                    const edge meshEdge0(meshPtsA, e0);
 
                     const auto iter = edgeValues.cfind(meshEdge0);
 
@@ -546,7 +548,7 @@ void Foam::syncTools::syncEdgeMap
                 }
                 {
                     const edge& e1 = edgesB[twoEdges[1]];
-                    const edge meshEdge1(meshPtsB[e1[0]], meshPtsB[e1[1]]);
+                    const edge meshEdge1(meshPtsB, e1);
 
                     const auto iter = edgeValues.cfind(meshEdge1);
 
@@ -573,7 +575,7 @@ void Foam::syncTools::syncEdgeMap
                 if (half1Fnd.good())
                 {
                     const edge& e0 = edgesA[twoEdges[0]];
-                    const edge meshEdge0(meshPtsA[e0[0]], meshPtsA[e0[1]]);
+                    const edge meshEdge0(meshPtsA, e0);
 
                     combine
                     (
@@ -589,7 +591,7 @@ void Foam::syncTools::syncEdgeMap
                 if (half0Fnd.good())
                 {
                     const edge& e1 = edgesB[twoEdges[1]];
-                    const edge meshEdge1(meshPtsB[e1[0]], meshPtsB[e1[1]]);
+                    const edge meshEdge1(meshPtsB, e1);
 
                     combine
                     (
@@ -749,8 +751,8 @@ void Foam::syncTools::syncPointList
     {
         FatalErrorInFunction
             << "Number of values " << pointValues.size()
-            << " is not equal to the number of points in the mesh "
-            << mesh.nPoints() << abort(FatalError);
+            << " != number of points " << mesh.nPoints() << nl
+            << abort(FatalError);
     }
 
     mesh.globalData().syncPointData(pointValues, cop, top);
@@ -772,8 +774,8 @@ void Foam::syncTools::syncPointList
     {
         FatalErrorInFunction
             << "Number of values " << pointValues.size()
-            << " is not equal to the number of meshPoints "
-            << meshPoints.size() << abort(FatalError);
+            << " != number of meshPoints " << meshPoints.size() << nl
+            << abort(FatalError);
     }
     const globalMeshData& gd = mesh.globalData();
     const indirectPrimitivePatch& cpp = gd.coupledPatch();
@@ -829,8 +831,8 @@ void Foam::syncTools::syncEdgeList
     {
         FatalErrorInFunction
             << "Number of values " << edgeValues.size()
-            << " is not equal to the number of edges in the mesh "
-            << mesh.nEdges() << abort(FatalError);
+            << " != number of edges " << mesh.nEdges() << nl
+            << abort(FatalError);
     }
 
     const edgeList& edges = mesh.edges();
@@ -915,7 +917,7 @@ template<class T, class CombineOp, class TransformOp, class FlipOp>
 void Foam::syncTools::syncEdgeList
 (
     const polyMesh& mesh,
-    const labelList& meshEdges,
+    const labelUList& meshEdges,
     List<T>& edgeValues,
     const CombineOp& cop,
     const T& nullValue,
@@ -927,8 +929,8 @@ void Foam::syncTools::syncEdgeList
     {
         FatalErrorInFunction
             << "Number of values " << edgeValues.size()
-            << " is not equal to the number of meshEdges "
-            << meshEdges.size() << abort(FatalError);
+            << " != number of meshEdges " << meshEdges.size() << nl
+            << abort(FatalError);
     }
     const edgeList& edges = mesh.edges();
     const globalMeshData& gd = mesh.globalData();
@@ -946,7 +948,7 @@ void Foam::syncTools::syncEdgeList
         const auto iter = mpm.cfind(meshEdgei);
         if (iter.good())
         {
-            const label cppEdgei = iter();
+            const label cppEdgei = iter.val();
             const edge& cppE = cppEdges[cppEdgei];
             const edge& meshE = edges[meshEdgei];
 
@@ -992,7 +994,7 @@ void Foam::syncTools::syncEdgeList
         const auto iter = mpm.cfind(meshEdgei);
         if (iter.good())
         {
-            label cppEdgei = iter();
+            label cppEdgei = iter.val();
             const edge& cppE = cppEdges[cppEdgei];
             const edge& meshE = edges[meshEdgei];
 
@@ -1029,14 +1031,13 @@ void Foam::syncTools::syncBoundaryFaceList
     {
         FatalErrorInFunction
             << "Number of values " << faceValues.size()
-            << " is not equal to the number of boundary faces in the mesh "
-            << mesh.nBoundaryFaces() << nl
+            << " != number of boundary faces " << mesh.nBoundaryFaces() << nl
             << abort(FatalError);
     }
 
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
-    if (parRun)
+    if (parRun && UPstream::parRun())
     {
         // Avoid mesh.globalData() - possible race condition
 
@@ -1104,7 +1105,7 @@ void Foam::syncTools::syncBoundaryFaceList
             }
 
             // Wait for all comms to finish
-            Pstream::waitRequests(startRequest);
+            UPstream::waitRequests(startRequest);
 
             // Combine with existing data
             for (const polyPatch& pp : patches)
@@ -1274,20 +1275,20 @@ void Foam::syncTools::syncFaceList
     {
         FatalErrorInFunction
             << "Number of values " << faceValues.size()
-            << " is not equal to the number of "
+            << " != number of "
             << (isBoundaryOnly ? "boundary" : "mesh") << " faces "
-            << ((mesh.nFaces() - boundaryOffset)) << nl
+            << (mesh.nFaces() - boundaryOffset) << nl
             << abort(FatalError);
     }
 
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
-    if (parRun)
+    if (parRun && UPstream::parRun())
     {
         const label startRequest = UPstream::nRequests();
 
         // Receive buffers
-        PtrList<PackedList<Width>> recvInfos(patches.size());
+        PtrList<PackedList<Width>> recvBufs(patches.size());
 
         // Set up reads
         for (const polyPatch& pp : patches)
@@ -1298,23 +1299,21 @@ void Foam::syncTools::syncFaceList
             {
                 const auto& procPatch = *ppp;
                 const label patchi = pp.index();
-                const label patchSize = pp.size();
 
-                recvInfos.set(patchi, new PackedList<Width>(patchSize));
-                PackedList<Width>& recvInfo = recvInfos[patchi];
+                auto& recvbuf = recvBufs.emplace_set(patchi, pp.size());
 
                 UIPstream::read
                 (
                     UPstream::commsTypes::nonBlocking,
                     procPatch.neighbProcNo(),
-                    recvInfo.data_bytes(),
-                    recvInfo.size_bytes()
+                    recvbuf.data_bytes(),
+                    recvbuf.size_bytes()
                 );
             }
         }
 
         // Send buffers
-        PtrList<PackedList<Width>> sendInfos(patches.size());
+        PtrList<PackedList<Width>> sendBufs(patches.size());
 
         // Set up writes
         for (const polyPatch& pp : patches)
@@ -1326,24 +1325,16 @@ void Foam::syncTools::syncFaceList
                 const auto& procPatch = *ppp;
                 const label patchi = pp.index();
 
-                const labelRange range
-                (
-                    pp.start()-boundaryOffset,
-                    pp.size()
-                );
-                sendInfos.set
-                (
-                    patchi,
-                    new PackedList<Width>(faceValues, range)
-                );
-                PackedList<Width>& sendInfo = sendInfos[patchi];
+                const labelRange range(pp.start()-boundaryOffset, pp.size());
+
+                auto& sendbuf = sendBufs.emplace_set(patchi, faceValues, range);
 
                 UOPstream::write
                 (
                     UPstream::commsTypes::nonBlocking,
                     procPatch.neighbProcNo(),
-                    sendInfo.cdata_bytes(),
-                    sendInfo.size_bytes()
+                    sendbuf.cdata_bytes(),
+                    sendbuf.size_bytes()
                 );
             }
         }
@@ -1361,13 +1352,13 @@ void Foam::syncTools::syncFaceList
                 const label patchi = pp.index();
                 const label patchSize = pp.size();
 
-                const PackedList<Width>& recvInfo = recvInfos[patchi];
+                const auto& recvbuf = recvBufs[patchi];
 
                 // Combine (bitwise)
                 label bFacei = pp.start()-boundaryOffset;
                 for (label i = 0; i < patchSize; ++i)
                 {
-                    unsigned int recvVal = recvInfo[i];
+                    unsigned int recvVal = recvbuf[i];
                     unsigned int faceVal = faceValues[bFacei];
 
                     cop(faceVal, recvVal);
@@ -1420,15 +1411,16 @@ void Foam::syncTools::swapBoundaryCellList
 (
     const polyMesh& mesh,
     const UList<T>& cellData,
-    List<T>& neighbourCellData
+    List<T>& neighbourCellData,
+    const bool parRun
 )
 {
     if (cellData.size() != mesh.nCells())
     {
         FatalErrorInFunction
             << "Number of cell values " << cellData.size()
-            << " is not equal to the number of cells in the mesh "
-            << mesh.nCells() << abort(FatalError);
+            << " != number of cells " << mesh.nCells() << nl
+            << abort(FatalError);
     }
 
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
@@ -1437,15 +1429,18 @@ void Foam::syncTools::swapBoundaryCellList
 
     for (const polyPatch& pp : patches)
     {
-        label bFacei = pp.offset();
+        const auto& faceCells = pp.faceCells();
 
-        for (const label celli : pp.faceCells())
-        {
-            neighbourCellData[bFacei] = cellData[celli];
-            ++bFacei;
-        }
+        // ie, boundarySlice() = patchInternalList()
+        SubList<T>
+        (
+            neighbourCellData,
+            faceCells.size(),
+            pp.offset()
+        ) = UIndirectList<T>(cellData, faceCells);
     }
-    syncTools::swapBoundaryFaceList(mesh, neighbourCellData);
+
+    syncTools::swapBoundaryFaceList(mesh, neighbourCellData, parRun);
 }
 
 
@@ -1479,10 +1474,17 @@ template<unsigned Width>
 void Foam::syncTools::swapFaceList
 (
     const polyMesh& mesh,
-    PackedList<Width>& faceValues
+    PackedList<Width>& faceValues,
+    const bool parRun
 )
 {
-    syncFaceList(mesh, faceValues, eqOp<unsigned int>());
+    syncFaceList
+    (
+        mesh,
+        faceValues,
+        eqOp<unsigned int>(),
+        parRun
+    );
 }
 
 
@@ -1490,10 +1492,17 @@ template<unsigned Width>
 void Foam::syncTools::swapBoundaryFaceList
 (
     const polyMesh& mesh,
-    PackedList<Width>& faceValues
+    PackedList<Width>& faceValues,
+    const bool parRun
 )
 {
-    syncBoundaryFaceList(mesh, faceValues, eqOp<unsigned int>());
+    syncBoundaryFaceList
+    (
+        mesh,
+        faceValues,
+        eqOp<unsigned int>(),
+        parRun
+    );
 }
 
 
@@ -1510,8 +1519,8 @@ void Foam::syncTools::syncPointList
     {
         FatalErrorInFunction
             << "Number of values " << pointValues.size()
-            << " is not equal to the number of points in the mesh "
-            << mesh.nPoints() << abort(FatalError);
+            << " != number of points " << mesh.nPoints() << nl
+            << abort(FatalError);
     }
 
     const globalMeshData& gd = mesh.globalData();
@@ -1553,8 +1562,8 @@ void Foam::syncTools::syncEdgeList
     {
         FatalErrorInFunction
             << "Number of values " << edgeValues.size()
-            << " is not equal to the number of edges in the mesh "
-            << mesh.nEdges() << abort(FatalError);
+            << " != number of edges " << mesh.nEdges() << nl
+            << abort(FatalError);
     }
 
     const globalMeshData& gd = mesh.globalData();

src/OpenFOAM/meshes/polyMesh/zones/ZoneMesh/ZoneMesh.C

@@ -96,7 +96,7 @@ void Foam::ZoneMesh<ZoneType, MeshType>::calcZoneMap() const
                 {
                     map.insert(id, zonei);
                 }
-                else if (fnd() != zonei)
+                else if (fnd.val() != zonei)
                 {
                     // Multiple zones for same id
 
@@ -112,7 +112,7 @@ void Foam::ZoneMesh<ZoneType, MeshType>::calcZoneMap() const
                               : static_cast<const labelList&>(zn)
                             )
                             {
-                                maxIndex = max(maxIndex, id);
+                                maxIndex = Foam::max(maxIndex, id);
                             }
                         }
                         additionalMapPtr_.reset(new labelListList(maxIndex+1));
@@ -135,7 +135,7 @@ void Foam::ZoneMesh<ZoneType, MeshType>::calcZoneMap() const
 
                 if (zones.size())
                 {
-                    stableSort(zones);
+                    Foam::stableSort(zones);
                     const label zonei = map[id];
                     const label index = findLower(zones, zonei);
                     if (index == -1)
@@ -426,11 +426,11 @@ Foam::label Foam::ZoneMesh<ZoneType, MeshType>::whichZones
 ) const
 {
     zones.clear();
-    const auto fnd = zoneMap().find(objectIndex);
+    const auto fnd = zoneMap().cfind(objectIndex);
     if (fnd)
     {
         // Add main element
-        zones.push_back(fnd());
+        zones.push_back(fnd.val());
         if (additionalMapPtr_)
         {
             const auto& additionalMap = *additionalMapPtr_;
@@ -967,7 +967,7 @@ bool Foam::ZoneMesh<ZoneType, MeshType>::checkParallelSync
     const bool report
 ) const
 {
-    if (!Pstream::parRun())
+    if (!UPstream::parRun())
     {
         return false;
     }
@@ -1025,7 +1025,7 @@ bool Foam::ZoneMesh<ZoneType, MeshType>::checkParallelSync
             {
                 hasError = true;
 
-                if (debug || (report && Pstream::master()))
+                if (debug || (report && UPstream::master()))
                 {
                     Info<< " ***Zone " << zn.name()
                         << " of type " << zn.type()

src/OpenFOAM/meshes/polyMesh/zones/ZoneMesh/ZoneMesh.H

@@ -194,9 +194,8 @@ public:
         //  If object does not belong to any zones, return -1
         label whichZone(const label objectIndex) const;
 
-        //- Given a global object index, return (in argument) the zones it is
-        //  in. Returns number of zones (0 if object does not belong to any
-        //  zones)
+        //- Given a global object index, return (in argument) its zones.
+        //  Returns number of zones (0 if object does not belong to any zones)
         label whichZones
         (
             const label objectIndex,

src/atmosphericModels/Make/files

@@ -33,7 +33,15 @@ fvOptions/atmPlantCanopyTurbSource/atmPlantCanopyTurbSource.C
 fvOptions/atmPlantCanopyUSource/atmPlantCanopyUSource.C
 fvOptions/atmPlantCanopyTSource/atmPlantCanopyTSource.C
 fvOptions/atmNutSource/atmNutSource.C
+fvOptions/treeTurbulence/treeTurbulence.C
 
+etht = fvOptions/evapotranspirationHeatTransfer
+ethtModels = $(etht)/evapotranspirationHeatTransferModels
+$(etht)/evapotranspirationHeatTransfer.C
+$(ethtModels)/evapotranspirationHeatTransferModel/evapotranspirationHeatTransferModel.C
+$(ethtModels)/evapotranspirationHeatTransferModel/evapotranspirationHeatTransferModelNew.C
+$(ethtModels)/tree/tree.C
+$(ethtModels)/grass/grass.C
 
 /* functionObjects */
 functionObjects/ObukhovLength/ObukhovLength.C

src/atmosphericModels/Make/options

@@ -8,6 +8,7 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
@@ -18,6 +19,7 @@ LIB_LIBS = \
     -lfvOptions \
     -lmeshTools \
     -lsurfMesh \
+    -lradiationModels \
     -lfluidThermophysicalModels \
     -lsolidThermo \
     -lturbulenceModels \

src/atmosphericModels/fvOptions/evapotranspirationHeatTransfer/evapotranspirationHeatTransfer.C

@@ -0,0 +1,172 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "evapotranspirationHeatTransfer.H"
+#include "evapotranspirationHeatTransferModel.H"
+#include "fvMesh.H"
+#include "fvMatrix.H"
+#include "basicThermo.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+    defineTypeNameAndDebug(evapotranspirationHeatTransfer, 0);
+    addToRunTimeSelectionTable
+    (
+        option,
+        evapotranspirationHeatTransfer,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fv::evapotranspirationHeatTransfer::evapotranspirationHeatTransfer
+(
+    const word& sourceName,
+    const word& modelType,
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    fv::cellSetOption(sourceName, modelType, dict, mesh),
+    ethtModelPtr_()
+{
+    // Set the field name to that of the energy
+    // field from which the temperature is obtained
+
+    const auto& thermo = mesh_.lookupObject<basicThermo>(basicThermo::dictName);
+
+    fieldNames_.resize(1, thermo.he().name());
+
+    fv::option::resetApplied();
+
+    Info<< "    Applying evapotranspirationHeatTransfer to: " << fieldNames_[0]
+        << endl;
+
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::fv::evapotranspirationHeatTransfer::~evapotranspirationHeatTransfer()
+{}  // evapotranspirationHeatTransferModel was forward declared
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::fv::evapotranspirationHeatTransfer::addSup
+(
+    fvScalarMatrix& eqn,
+    const label fieldi
+)
+{
+    if (this->V() < VSMALL)
+    {
+        return;
+    }
+
+    const scalar V = this->V();
+    const scalarField& Vcells = mesh_.V();
+
+    const volScalarField& he = eqn.psi();
+    scalarField& heSource = eqn.source();
+
+    // Calculate evapotranspiration heat transfer rate per volume [J/s/m^3]
+    const scalarField Q(ethtModelPtr_->Q(cells_)/V);
+
+    if (he.dimensions() == dimEnergy/dimMass)
+    {
+        forAll(cells_, i)
+        {
+            const label celli = cells_[i];
+
+            heSource[celli] += Q[i]*Vcells[celli];
+        }
+    }
+    else if (he.dimensions() == dimTemperature)
+    {
+        const auto& thermo =
+            mesh_.lookupObject<basicThermo>(basicThermo::dictName);
+
+        // Calculate density*heat capacity at constant pressure/volume
+        const volScalarField rhoCpv(thermo.rho()*thermo.Cpv());
+
+        // heSource [K m^3/s] = [J/s/m^3] * m^3 / [kg/m^3] / [J/kg/K]
+        forAll(cells_, i)
+        {
+            const label celli = cells_[i];
+
+            heSource[celli] += Q[i]*Vcells[celli]*rhoCpv[i];
+        }
+    }
+}
+
+
+void Foam::fv::evapotranspirationHeatTransfer::addSup
+(
+    const volScalarField& rho,
+    fvScalarMatrix& eqn,
+    const label fieldi
+)
+{
+    addSup(eqn, fieldi);
+}
+
+
+bool Foam::fv::evapotranspirationHeatTransfer::read(const dictionary& dict)
+{
+    if (!fv::cellSetOption::read(dict))
+    {
+        return false;
+    }
+
+    if (selectionMode_ != selectionModeType::smCellZone)
+    {
+        FatalIOErrorInFunction(dict)
+            << "evapotranspirationHeatTransfer requires "
+            << selectionModeTypeNames_[selectionModeType::smCellZone]
+            << exit(FatalIOError);
+    }
+
+    ethtModelPtr_.reset
+    (
+        evapotranspirationHeatTransferModel::New(dict, mesh_)
+    );
+
+    return true;
+}
+
+
+// ************************************************************************* //

src/atmosphericModels/fvOptions/evapotranspirationHeatTransfer/evapotranspirationHeatTransfer.H

@@ -0,0 +1,183 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::fv::evapotranspirationHeatTransfer
+
+Description
+    Applies sources on temperature (\c T - incompressible)
+    or energy (\c h/e - compressible) equation to incorporate
+    evapotranspiration heat-transfer effects from the specified
+    plant canopy. Heat transfer is usually calculated based on
+    empirical relations between plants and solar radiation.
+
+Usage
+    Example by using \c constant/fvOptions:
+    \verbatim
+    evapotranspirationHeatTransfer1
+    {
+        // Mandatory entries
+        type                evapotranspirationHeatTransfer;
+        model               <word>;
+
+        // Model-specific entries
+        ...
+
+        // Inherited entries
+        ...
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property | Description                         | Type | Reqd | Deflt
+      type     | Type name: evapotranspirationHeatTransfer  | word | yes  | -
+      model    | Name of heat-transfer model        | word  | yes  | -
+    \endtable
+
+    Options for the \c model entry:
+    \verbatim
+      tree   | Regression-based heat-transfer model for trees
+      grass  | Heat-transfer model for grass
+    \endverbatim
+
+    The inherited entries are elaborated in:
+      - \link fvOption.H \endlink
+      - \link cellSetOption.H \endlink
+
+Note
+  - Evapotranspiration mass transfer is not modelled.
+
+SourceFiles
+    evapotranspirationHeatTransfer.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_fv_evapotranspirationHeatTransfer_H
+#define Foam_fv_evapotranspirationHeatTransfer_H
+
+#include "cellSetOption.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+class evapotranspirationHeatTransferModel;
+
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+                Class evapotranspirationHeatTransfer Declaration
+\*---------------------------------------------------------------------------*/
+
+class evapotranspirationHeatTransfer
+:
+    public fv::cellSetOption
+{
+    // Private Data
+
+        //- Runtime-selectable evapotranspiration heat transfer model
+        autoPtr<evapotranspirationHeatTransferModel> ethtModelPtr_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("evapotranspirationHeatTransfer");
+
+
+    // Constructors
+
+        //- Construct from explicit source name and mesh
+        evapotranspirationHeatTransfer
+        (
+            const word& sourceName,
+            const word& modelType,
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+        //- No copy construct
+        evapotranspirationHeatTransfer(const evapotranspirationHeatTransfer&) =
+            delete;
+
+        //- No copy assignment
+        void operator=(const evapotranspirationHeatTransfer&) = delete;
+
+
+    //- Destructor
+    virtual ~evapotranspirationHeatTransfer();
+
+
+    // Member Functions
+
+        //- Add explicit contribution to energy equation
+        //- for incompressible flow computations
+        virtual void addSup
+        (
+            fvScalarMatrix& eqn,
+            const label fieldi
+        );
+
+        //- Add explicit contribution to energy equation
+        //- for compressible flow computations
+        virtual void addSup
+        (
+            const volScalarField& rho,
+            fvScalarMatrix& eqn,
+            const label fieldi
+        );
+
+        //- Add explicit contribution to energy equation
+        //- for multiphase flow computations
+        virtual void addSup
+        (
+            const volScalarField& alpha,
+            const volScalarField& rho,
+            fvScalarMatrix& eqn,
+            const label fieldi
+        )
+        {
+            NotImplemented;
+        }
+
+        //- Read dictionary
+        virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/atmosphericModels/fvOptions/evapotranspirationHeatTransfer/evapotranspirationHeatTransferModels/evapotranspirationHeatTransferModel/evapotranspirationHeatTransferModel.C

@@ -0,0 +1,196 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "evapotranspirationHeatTransferModel.H"
+#include "radiationModel.H"
+#include "solarLoadBase.H"
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(evapotranspirationHeatTransferModel, 0);
+    defineRunTimeSelectionTable
+    (
+        evapotranspirationHeatTransferModel,
+        dictionary
+    );
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+Foam::volScalarField& Foam::evapotranspirationHeatTransferModel::getOrReadField
+(
+    const word& fieldName
+) const
+{
+    auto* ptr = mesh_.getObjectPtr<volScalarField>(fieldName);
+
+    if (!ptr)
+    {
+        ptr = new volScalarField
+        (
+            IOobject
+            (
+                fieldName,
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::MUST_READ,
+                IOobject::AUTO_WRITE
+            ),
+            mesh_
+        );
+        mesh_.objectRegistry::store(ptr);
+    }
+
+    return *ptr;
+}
+
+
+Foam::scalar Foam::evapotranspirationHeatTransferModel::q() const
+{
+    // Retrieve solar-load information
+    const auto& base =
+        mesh().lookupObject<radiation::solarLoadBase>("solarLoadBase");
+    const solarCalculator& solar = base.solarCalculatorRef();
+
+    // Retrieve internal & boundary faces directly hit by solar rays
+    const faceShading& shade = base.faceShadingRef();
+    const labelList& hitFacesId = shade.rayStartFaces();
+
+    // Retrieve face zone information
+    const faceZone& zone = mesh_.faceZones()[faceZoneId_];
+    const vectorField& faceNormals = zone().faceNormals();
+    const scalarField& faceAreas = zone().magFaceAreas();
+
+    // Retrieve direct solar load [W/m^2]
+    const vector directSolarRad(solar.directSolarRad()*solar.direction());
+
+    // Identify zone faces within mesh
+    bitSet isZoneFace(mesh_.nFaces());
+    isZoneFace.set(zone);
+
+    // Calculate area-averaged incident solar load
+    // Assuming hit faces are updated by solarLoad
+    scalar q = 0;
+    scalar totalFaceArea = 0;
+
+    for (const label facei : hitFacesId)
+    {
+        if (isZoneFace[facei])
+        {
+            const label localFacei = zone.whichFace(facei);
+
+            const vector& faceNormal = faceNormals[localFacei];
+            const scalar faceArea = faceAreas[localFacei];
+
+            const scalar qIncident = directSolarRad & faceNormal;
+
+            q += qIncident*faceArea;
+            totalFaceArea += faceArea;
+        }
+    }
+    reduce(q, sumOp<scalar>());
+    reduce(totalFaceArea, sumOp<scalar>());
+    q /= totalFaceArea;
+
+
+    // Sum diffusive solar loads
+    switch(solar.sunLoadModel())
+    {
+        case solarCalculator::mSunLoadFairWeatherConditions:
+        case solarCalculator::mSunLoadTheoreticalMaximum:
+        {
+            // Calculate diffusive radiance
+            // contribution from sky and ground
+            tmp<scalarField> tdiffuseSolarRad =
+                solar.diffuseSolarRad(faceNormals);
+            const scalarField& diffuseSolarRad = tdiffuseSolarRad.cref();
+
+            // Calculate area-averaged diffusive solar load
+            scalar meanDiffuseSolarRad = 0;
+            forAll(faceAreas, i)
+            {
+                meanDiffuseSolarRad += diffuseSolarRad[i]*faceAreas[i];
+            }
+            meanDiffuseSolarRad /= totalFaceArea;
+
+            q += meanDiffuseSolarRad;
+        }
+        break;
+
+        case solarCalculator::mSunLoadConstant:
+        case solarCalculator::mSunLoadTimeDependent:
+        {
+            q += solar.diffuseSolarRad();
+        }
+        break;
+    }
+
+    return q;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::evapotranspirationHeatTransferModel::evapotranspirationHeatTransferModel
+(
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    mesh_(mesh)
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+bool Foam::evapotranspirationHeatTransferModel::read(const dictionary& dict)
+{
+    faceZoneId_ = mesh_.faceZones().findZoneID(dict.get<word>("faceZone"));
+
+    if (faceZoneId_ < 0)
+    {
+        const word faceZoneName(mesh_.faceZones()[faceZoneId_].name());
+
+        FatalIOErrorInFunction(dict)
+            << type() << ' ' << typeName << ": "
+            << "    No matching face zone: " << faceZoneName  << nl
+            << "    Known face zones: "
+            << flatOutput(mesh_.faceZones().names()) << nl
+            << exit(FatalIOError);
+
+        return false;
+    }
+
+    return true;
+}
+
+
+// ************************************************************************* //

src/atmosphericModels/fvOptions/evapotranspirationHeatTransfer/evapotranspirationHeatTransferModels/evapotranspirationHeatTransferModel/evapotranspirationHeatTransferModel.H

@@ -0,0 +1,180 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::evapotranspirationHeatTransferModel
+
+Description
+    Base class for evapotranspiration models
+    to handle general evapotranspiration characteristics.
+
+SourceFiles
+    evapotranspirationHeatTransferModel.C
+    evapotranspirationHeatTransferModelNew.C
+    evapotranspirationHeatTransferModelTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_evapotranspirationHeatTransferModel_H
+#define Foam_evapotranspirationHeatTransferModel_H
+
+#include "dictionary.H"
+#include "HashSet.H"
+#include "volFields.H"
+#include "runTimeSelectionTables.H"
+#include "OFstream.H"
+#include "coordinateSystem.H"
+#include "writeFile.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+class fvMesh;
+
+/*---------------------------------------------------------------------------*\
+            Class evapotranspirationHeatTransferModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class evapotranspirationHeatTransferModel
+{
+    // Private Data
+
+        //- Reference to the mesh
+        const fvMesh& mesh_;
+
+        //- Identifier of specified face zone
+        label faceZoneId_;
+
+
+protected:
+
+    // Protected Member Functions
+
+        //- Return requested field from the object registry
+        //- or read+register the field to the object registry
+        volScalarField& getOrReadField(const word& fieldName) const;
+
+        //- Area-averaged heat flux due to short-wave
+        //- solar rays on face zone [W/m^2]
+        //  Short-wave: direct and diffusive solar rays
+        scalar q() const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("evapotranspirationHeatTransferModel");
+
+
+    // Declare runtime constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            evapotranspirationHeatTransferModel,
+            dictionary,
+            (
+                const dictionary& dict,
+                const fvMesh& mesh
+            ),
+            (dict, mesh)
+        );
+
+
+    // Selectors
+
+        //- Return a reference to the selected evapotranspiration model
+        static autoPtr<evapotranspirationHeatTransferModel> New
+        (
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+
+    // Constructors
+
+        //- Construct from components
+        evapotranspirationHeatTransferModel
+        (
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+        //- No copy construct
+        evapotranspirationHeatTransferModel
+        (
+            const evapotranspirationHeatTransferModel&
+        ) = delete;
+
+        //- No copy assignment
+        void operator=(const evapotranspirationHeatTransferModel&) = delete;
+
+
+    //- Destructor
+    virtual ~evapotranspirationHeatTransferModel() = default;
+
+
+    // Member Functions
+
+    // Access
+
+        //- Return const reference to the mesh
+        const fvMesh& mesh() const noexcept
+        {
+            return mesh_;
+        }
+
+        //- Return the face-zone identifier
+        label faceZoneId() const noexcept
+        {
+            return faceZoneId_;
+        }
+
+
+    // Evaluation
+
+        //- Return heat-transfer rate for speficied cells [J/s]
+        virtual tmp<scalarField> Q(const labelList& cells) const = 0;
+
+
+    // I-O
+
+        //- Read the dictionary
+        virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/atmosphericModels/fvOptions/evapotranspirationHeatTransfer/evapotranspirationHeatTransferModels/evapotranspirationHeatTransferModel/evapotranspirationHeatTransferModelNew.C

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

src/atmosphericModels/fvOptions/evapotranspirationHeatTransfer/evapotranspirationHeatTransferModels/grass/grass.C

@@ -0,0 +1,418 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "grass.H"
+#include "basicThermo.H"
+#include "fluidThermo.H"
+#include "turbulentTransportModel.H"
+#include "turbulentFluidThermoModel.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace evapotranspirationHeatTransferModels
+{
+    defineTypeNameAndDebug(grass, 0);
+    addToRunTimeSelectionTable
+    (
+        evapotranspirationHeatTransferModel,
+        grass,
+        dictionary
+    );
+}
+}
+
+
+const Foam::Enum
+<
+    Foam::evapotranspirationHeatTransferModels::grass::soilHeatFluxType
+>
+Foam::evapotranspirationHeatTransferModels::grass::soilHeatFluxTypeNames
+({
+    {
+        soilHeatFluxType::PROPORTIONAL_TO_SOLAR_RADIATION,
+        "proportionalToSolarRadiation"
+    },
+    { soilHeatFluxType::BOUNDARY, "boundary" }
+});
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::tmp<Foam::scalarField>
+Foam::evapotranspirationHeatTransferModels::grass::E(const labelList& cells)
+const
+{
+    const scalar q = this->q();
+    const scalar G = this->G(cells);
+    const scalar Delta = this->Delta();
+    const scalar D = this->D();
+    const scalar ra = this->ra();
+    const scalar rs = this->rs();
+    const scalar gamma = this->gamma();
+
+    // (BSG:Eq. 11)
+    const scalar E =
+        (Delta*(q - G) + rho_*Cp_*D/ra)/(Delta + gamma*(1 + rs/ra));
+
+    return tmp<scalarField>::New(cells.size(), E);
+}
+
+
+Foam::scalar Foam::evapotranspirationHeatTransferModels::grass::Delta() const
+{
+    // (BSG:Eq. 15), (APR:Eq. 13) - note 0.6108 -> 610.8 due to kPa -> Pa
+    const scalar Ta = Tref_ + 237.3;
+
+    return 4098*(610.8*exp(17.27*Tref_/Ta))/sqr(Ta);
+}
+
+
+Foam::scalar Foam::evapotranspirationHeatTransferModels::grass::D() const
+{
+    // (BSG:Eq. 16)
+    const scalar RHref = RHptr_->value(mesh().time().timeOutputValue());
+
+    return (1 - RHref)*pSat();
+}
+
+
+Foam::scalar Foam::evapotranspirationHeatTransferModels::grass::pSat() const
+{
+    // (BSG:Eq. 17; Arden Buck equation)
+    const scalar p1 = (1.0007 + 3.46e-8*pAtm_)*611.21;
+
+    return p1*exp((18.678 - Tref_/234.5)*Tref_/(Tref_ + 257.14));
+}
+
+
+Foam::scalar Foam::evapotranspirationHeatTransferModels::grass::gamma() const
+{
+    // (APR:Eq. 8)
+    return Cp_*pAtm_/(epsilon_*lambda());
+}
+
+
+Foam::scalar Foam::evapotranspirationHeatTransferModels::grass::lambda() const
+{
+    // (BSG:Eq. 12)
+    const scalar T1 = Tref_ + 273.15;
+    const scalar T2 = Tref_ + 239.24;
+
+    return 1.91846e6*sqr(T1/T2);
+}
+
+
+Foam::scalar Foam::evapotranspirationHeatTransferModels::grass::ra() const
+{
+    // (APR:Eq. 4), (BSG:Eq. 14)
+    const scalar log1 = log((zRefU_ - d_*h_)/(zom_*h_));
+    const scalar log2 = log((zRefH_ - d_*h_)/(zoh_*h_));
+
+    return log1*log2/(sqr(kappa_)*uRef_);
+}
+
+
+Foam::scalar Foam::evapotranspirationHeatTransferModels::grass::rs() const
+{
+    // (BSG:Eq. 13), (APR:Eq. 5; reduced form in p. 4-5)
+    return ri_/(scalar(12)*h_);
+}
+
+
+Foam::scalar Foam::evapotranspirationHeatTransferModels::grass::S
+(
+    const labelList& cells
+) const
+{
+    // Mark fvOption cells within mesh
+    bitSet isZoneCell(mesh().nCells());
+    isZoneCell.set(cells);
+
+    scalar S = 0;
+
+    // Select cells next to specified patches
+    // Sum patch area that is covered by fvOption cells
+    for (const label patchi : patchSet_)
+    {
+        const scalarField& s = mesh().magSf().boundaryField()[patchi];
+
+        const polyPatch& pp = mesh().boundaryMesh()[patchi];
+        const labelList& faceCells = pp.faceCells();
+
+        forAll(faceCells, i)
+        {
+            const bool isCovered = isZoneCell[faceCells[i]];
+
+            if (isCovered)
+            {
+                S += s[i];
+            }
+        }
+    }
+    reduce(S, sumOp<scalar>());
+
+    if (mag(S) < SMALL)
+    {
+        FatalErrorInFunction
+            << "Area coverage of grass cannot be zero. "
+            << "Check whether cellZone has any boundary faces."
+            << exit(FatalError);
+    }
+
+    return S;
+}
+
+
+Foam::scalar Foam::evapotranspirationHeatTransferModels::grass::G
+(
+    const labelList& cells
+) const
+{
+    switch (soilHeatFluxMethod_)
+    {
+        case soilHeatFluxType::PROPORTIONAL_TO_SOLAR_RADIATION:
+        {
+            return Csoil_*q();
+        }
+
+        case soilHeatFluxType::BOUNDARY:
+        {
+            // Retrieve heat flux through patches
+            tmp<FieldField<Field, scalar>> tqBf = this->qBf();
+            const FieldField<Field, scalar>& qBf = tqBf.cref();
+
+            // Mark fvOption cells within mesh
+            bitSet isZoneCell(mesh().nCells());
+            isZoneCell.set(cells);
+
+            scalar G = 0;
+
+            for (const label patchi : patchSet_)
+            {
+                const scalarField& s = mesh().magSf().boundaryField()[patchi];
+
+                const scalarField& qfld = qBf[patchi];
+
+                const polyPatch& pp = mesh().boundaryMesh()[patchi];
+                const labelList& faceCells = pp.faceCells();
+
+                forAll(faceCells, i)
+                {
+                    const bool isCovered = isZoneCell[faceCells[i]];
+
+                    if (isCovered)
+                    {
+                        G += qfld[i]*s[i];
+                    }
+                }
+            }
+            reduce(G, sumOp<scalar>());
+
+            return G/area_;
+        }
+    }
+
+    return -1;
+}
+
+
+Foam::tmp<Foam::FieldField<Foam::Field, Foam::scalar>>
+Foam::evapotranspirationHeatTransferModels::grass::qBf() const
+{
+    const auto& T = mesh().lookupObject<volScalarField>(TName_);
+    const volScalarField::Boundary& Tbf = T.boundaryField();
+
+    auto tq = tmp<FieldField<Field, scalar>>::New(Tbf.size());
+    auto& q = tq.ref();
+
+    forAll(q, patchi)
+    {
+        q.set(patchi, new Field<scalar>(Tbf[patchi].size(), Zero));
+    }
+
+    typedef compressible::turbulenceModel cmpTurbModel;
+
+    if (mesh().foundObject<cmpTurbModel>(cmpTurbModel::propertiesName))
+    {
+        const auto& turb =
+            mesh().lookupObject<cmpTurbModel>(cmpTurbModel::propertiesName);
+
+        // Note: calling he(p,T) instead of he()
+        const volScalarField he(turb.transport().he(turb.transport().p(), T));
+        const volScalarField::Boundary& hebf = he.boundaryField();
+
+        const volScalarField alphaEff(turb.alphaEff());
+        const volScalarField::Boundary& alphaEffbf = alphaEff.boundaryField();
+
+        for (const label patchi : patchSet_)
+        {
+            q[patchi] = alphaEffbf[patchi]*hebf[patchi].snGrad();
+        }
+    }
+    else if (mesh().foundObject<fluidThermo>(fluidThermo::dictName))
+    {
+        const auto& thermo =
+            mesh().lookupObject<fluidThermo>(fluidThermo::dictName);
+
+        // Note: calling he(p,T) instead of he()
+        const volScalarField he(thermo.he(thermo.p(), T));
+        const volScalarField::Boundary& hebf = he.boundaryField();
+
+        const volScalarField& alpha(thermo.alpha());
+        const volScalarField::Boundary& alphabf = alpha.boundaryField();
+
+        for (const label patchi : patchSet_)
+        {
+            q[patchi] = alphabf[patchi]*hebf[patchi].snGrad();
+        }
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Unable to find a valid thermo model to evaluate q. " << nl
+            << "Database contents are: " << mesh().objectRegistry::sortedToc()
+            << exit(FatalError);
+    }
+
+    // No radiative heat flux contribution is present
+
+    return tq;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::evapotranspirationHeatTransferModels::grass::grass
+(
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    evapotranspirationHeatTransferModel(dict, mesh),
+    soilHeatFluxMethod_
+    (
+        soilHeatFluxTypeNames.getOrDefault
+        (
+            "soilHeatFluxMethod",
+            dict,
+            soilHeatFluxType::PROPORTIONAL_TO_SOLAR_RADIATION
+        )
+    ),
+    Tptr_(nullptr),
+    RHptr_(nullptr),
+    area_(-1),
+    Csoil_(),
+    rho_(),
+    Cp_(),
+    epsilon_(),
+    h_(),
+    kappa_(),
+    uRef_(),
+    zRefU_(),
+    zRefH_(),
+    zom_(),
+    zoh_(),
+    d_(),
+    ri_(),
+    pAtm_(),
+    Tref_(0),
+    TName_(),
+    patchSet_()
+{
+    Info<< "    Activating evapotranspiration heat transfer model: "
+        << typeName << endl;
+
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::scalarField>
+Foam::evapotranspirationHeatTransferModels::grass::Q
+(
+    const labelList& cells
+) const
+{
+    if (area_ == -1)
+    {
+        area_ = S(cells);
+    }
+
+    Tref_ = Tptr_->value(mesh().time().timeOutputValue());
+
+    return E(cells)*area_;
+}
+
+
+bool Foam::evapotranspirationHeatTransferModels::grass::read
+(
+    const dictionary& dict
+)
+{
+    if (!evapotranspirationHeatTransferModel::read(dict))
+    {
+        return false;
+    }
+
+    Tptr_.reset(Function1<scalar>::New("Tref", dict, &mesh()));
+    RHptr_.reset(Function1<scalar>::New("RHref", dict, &mesh()));
+
+    area_ = -1;
+
+    const auto range = scalarMinMax::ge(SMALL);
+
+    Csoil_ = dict.getCheckOrDefault<scalar>("Csoil", 0.1, range);
+    rho_ = dict.getCheckOrDefault<scalar>("rho", 1.225, range);
+    Cp_ = dict.getCheckOrDefault<scalar>("Cp", 1013.0, range);
+    epsilon_ = dict.getCheckOrDefault<scalar>("epsilon", 0.622, range);
+    h_ = dict.getCheckOrDefault<scalar>("h", 0.1, range);
+    kappa_ = dict.getCheckOrDefault<scalar>("kappa", 0.41, range);
+    uRef_ = dict.getCheckOrDefault<scalar>("uRef", 2, range);
+    zRefU_ = dict.getCheckOrDefault<scalar>("zRefU", 10, range);
+    zRefH_ = dict.getCheckOrDefault<scalar>("zRefH", 10, range);
+    zom_ = dict.getCheckOrDefault<scalar>("zom", 0.123, range);
+    zoh_ = dict.getCheckOrDefault<scalar>("zoh", 0.0123, range);
+    d_ = dict.getCheckOrDefault<scalar>("d", 2.0/3.0, range);
+    ri_ = dict.getCheckOrDefault<scalar>("ri", 100, range);
+    pAtm_ = dict.getCheckOrDefault<scalar>("pAtm", 101.325, range);
+    TName_ = dict.getOrDefault<word>("T", "T");
+
+    if (soilHeatFluxMethod_ == soilHeatFluxType::BOUNDARY)
+    {
+        patchSet_ =
+            mesh().boundaryMesh().patchSet(dict.get<wordRes>("patches"));
+    }
+
+    return true;
+}
+
+
+// ************************************************************************* //

src/atmosphericModels/fvOptions/evapotranspirationHeatTransfer/evapotranspirationHeatTransferModels/grass/grass.H

@@ -0,0 +1,342 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::evapotranspirationHeatTransferModels::grass
+
+Description
+    Applies sources on temperature (\c T - incompressible) or energy
+    (\c h/e - compressible) equation to incorporate evapotranspiration
+    heat-transfer effects from the specified grass canopy.
+    The model is based on Pemnan-Monteith Equation model.
+
+    Sources applied to either of the below, if exist:
+    \verbatim
+      T         | Temperature                                [K]
+      e         | Internal energy                            [m^2/s^2]
+      h         | Enthalphy                                  [m^2/s^2]
+    \endverbatim
+
+    Required fields:
+    \verbatim
+      T         | Temperature                                [K]
+      e         | Internal energy                            [m^2/s^2]
+      h         | Enthalphy                                  [m^2/s^2]
+      LAD       | Leaf area density                          [m^2/m^3]
+    \endverbatim
+
+    References:
+    \verbatim
+        Governing equations (tag:BSG):
+            Brozovsky, J., Simonsen, A., & Gaitani, N. (2021).
+            Validation of a CFD model for the evaluation of urban
+            microclimate at high latitudes: A case study in Trondheim, Norway.
+            Building and Environment, 205, 108175.
+            DOI:10.1016/j.buildenv.2021.108175
+
+        Governing equations (tag:APR):
+            Allen, R. G., Pereira, L. S., Raes, D., & Smith, M. (1998).
+            Crop evapotranspiration-Guidelines for computing crop
+            water requirements-FAO Irrigation and drainage paper 56.
+            Fao, Rome, 300(9), D05109.
+    \endverbatim
+
+Usage
+    Example by using \c constant/fvOptions:
+    \verbatim
+    evapotranspirationHeatTransfer1
+    {
+        // Inherited entries
+        ...
+
+        // Mandatory entries
+        model               grass;
+        Tref                <Function1<scalar>>;
+        RHref               <Function1<scalar>>;
+
+        // Conditional entries
+
+            // when soilHeatFluxMethod == boundary
+            patches             <wordRes>;
+
+        // Optional entries
+        soilHeatFluxMethod  <word>;
+        Csoil               <scalar>;
+        rho                 <scalar>;
+        Cp                  <scalar>;
+        epsilon             <scalar>;
+        h                   <scalar>;
+        kappa               <scalar>;
+        uRef                <scalar>;
+        zRefU               <scalar>;
+        zRefH               <scalar>;
+        zom                 <scalar>;
+        zoh                 <scalar>;
+        d                   <scalar>;
+        ri                  <scalar>;
+        pAtm                <scalar>;
+        T                   <word>;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property | Description                         | Type | Reqd | Deflt
+      model    | Model name: grass                   | word | yes  | -
+      Tref     | Reference weather station air temperature [Celsius] <!--
+               -->                     | Function1\<scalar\> | yes | -
+      RHref    | Reference weather station relative humidity [%]     <!--
+               -->                     | Function1\<scalar\> | yes | -
+      patches  | Names of ground patches           | wordRes | yes | -
+      soilHeatFluxMethod | Method to calculate soil heat flux - see below <!--
+               -->                   | word | no | -
+      Csoil    | Proportionality constant of soil heat-flux wrt <!--
+               --> solar radiation                   | scalar | no | 0.1
+      rho      | Mean air density at constant pressure [kg/m^3]      <!--
+               -->                                   | scalar | no | 1.225
+      Cp       | Specific heat at constant pressure [J/kg/C]         <!--
+               -->                                   | scalar | no | 1013.0
+      epsilon  | Molecular-weight ratio of water vapour/dry air [-]  <!--
+               -->                                  | scalar | no  | 0.622
+      h        | Height of grass layer [m]          | scalar | no  | 0.1
+      kappa    | Von Karman constant [-]            | scalar | no  | 0.41
+      uRef     | Reference velocity magnitude [m/s] | scalar | no  | 2.0
+      zRefU    | Height of wind speed measurements [m] | scalar | no | 10.0
+      zRefH    | Height of humidity measurements [m]   | scalar | no | 10.0
+      zom      | Roughness length governing momentum transfer coefficient <!--
+               -->                                  [-] | scalar | no | 0.123
+      zoh      | Roughness length governing transfer of heat and vapour <!--
+               --> coefficient [-]                      | scalar | no | 0.0123
+      d        | Zero plane displacement height coefficient [-]         <!--
+               -->                                      | scalar | no | 0.666
+      ri       | Bulk stomata resistance of leaf [s/m] | scalar | no | 100.0
+      pAtm     | Atmospheric pressure [kPa]            | scalar | no | 101.325
+      T        | Name of temperature field             | word   | no | T
+    \endtable
+
+    Options for the \c soilHeatFluxMethod entry:
+    \verbatim
+      proportionalToSolarRadiation | Estimate from solar load
+      boundary                     | Obtain soil heat flux from boundary
+    \endverbatim
+
+    The inherited entries are elaborated in:
+      - \link evapotranspirationHeatTransfer.H \endlink
+      - \link evapotranspirationHeatTransferModel.H \endlink
+
+SourceFiles
+    grass.C
+    grassTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_evapotranspirationHeatTransferModels_grass_H
+#define Foam_evapotranspirationHeatTransferModels_grass_H
+
+#include "evapotranspirationHeatTransferModel.H"
+#include "Function1.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace evapotranspirationHeatTransferModels
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class grass Declaration
+\*---------------------------------------------------------------------------*/
+
+class grass
+:
+    public evapotranspirationHeatTransferModel
+{
+    // Private Enumerations
+
+        //- Options for the soil heat flux
+        enum soilHeatFluxType : char
+        {
+            PROPORTIONAL_TO_SOLAR_RADIATION = 0, //!< "Estimate from solar load"
+            BOUNDARY,            //!< "Obtain soil heat flux from boundary"
+        };
+
+        //- Names for soilHeatFluxType
+        static const Enum<soilHeatFluxType> soilHeatFluxTypeNames;
+
+
+    // Private Data
+
+        //- Soil heat-flux calculation method
+        enum soilHeatFluxType soilHeatFluxMethod_;
+
+        //- Reference weather station air temperature [Celsius]
+        autoPtr<Function1<scalar>> Tptr_;
+
+        //- Reference weather station relative humidity [%]
+        autoPtr<Function1<scalar>> RHptr_;
+
+        //- Area coverage of grass [m^2]
+        mutable scalar area_;
+
+        //- Proportionality constant of soil heat-flux wrt solar radiation [-]
+        scalar Csoil_;
+
+        //- Mean air density at constant pressure [kg/m^3]
+        scalar rho_;
+
+        //- Specific heat at constant pressure [J/kg/C]
+        scalar Cp_;
+
+        //- Molecular-weight ratio of water vapour/dry air [-]
+        scalar epsilon_;
+
+        //- Height of grass layer [m]
+        scalar h_;
+
+        //- Von Karman constant [-]
+        scalar kappa_;
+
+        //- Reference velocity magnitude [m/s]
+        scalar uRef_;
+
+        //- Height of wind speed measurements [m]
+        scalar zRefU_;
+
+        //- Height of humidity measurements [m]
+        scalar zRefH_;
+
+        //- Roughness length governing momentum transfer coefficient [-]
+        scalar zom_;
+
+        //- Roughness length governing transfer of heat and vapour coeff [-]
+        scalar zoh_;
+
+        //- Zero plane displacement height coefficient [-]
+        scalar d_;
+
+        //- Bulk stomata resistance of leaf [s/m]
+        scalar ri_;
+
+        //- Atmospheric pressure [kPa]
+        scalar pAtm_;
+
+        //- Cached reference weather station air temperature [Celsius]
+        mutable scalar Tref_;
+
+        //- Name of temperature field
+        word TName_;
+
+        //- List of patches to calculate boundary heat flux
+        labelHashSet patchSet_;
+
+
+    // Private Member Functions
+
+        //- Return evapotranspiration heat flux through grass layer [W/m^2]
+        tmp<scalarField> E(const labelList& cells) const;
+
+        //- Return slope of the relation between
+        //- vapour pressure-temperature [Pa/K]
+        scalar Delta() const;
+
+        //- Return atmospheric vapour pressure deficit [Pa]
+        scalar D() const;
+
+        //- Return saturated vapour pressure over water [Pa]
+        scalar pSat() const;
+
+        //- Return psychrometric constant [kPa/K]
+        scalar gamma() const;
+
+        //- Return specific latent heat of vaporisation [MJ/kg]
+        scalar lambda() const;
+
+        //- Return bulk aerodynamic resistance [Pa]
+        scalar ra() const;
+
+        //- Return bulk surface resistance [Pa]
+        scalar rs() const;
+
+        //- Return area coverage of grass [m^2]
+        scalar S(const labelList& cells) const;
+
+        //- Return area-averaged heat flux through ground [W/m^2]
+        scalar G(const labelList& cells) const;
+
+        //- Return heat-flux boundary fields
+        tmp<FieldField<Field, scalar>> qBf() const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("grass");
+
+
+    // Constructors
+
+        //- Construct from components
+        grass
+        (
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+        //- No copy construct
+        grass(const grass&) = delete;
+
+        //- No copy assignment
+        void operator=(const grass&) = delete;
+
+
+    //- Destructor
+    virtual ~grass() = default;
+
+
+    // Member Functions
+
+    // Evaluation
+
+        //- Return heat-transfer rate for speficied cells [J/s]
+        virtual tmp<scalarField> Q(const labelList& cells) const;
+
+
+    // I-O
+
+        //- Read the dictionary
+        virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace evapotranspirationHeatTransferModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/atmosphericModels/fvOptions/evapotranspirationHeatTransfer/evapotranspirationHeatTransferModels/tree/tree.C

@@ -0,0 +1,139 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "tree.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace evapotranspirationHeatTransferModels
+{
+    defineTypeNameAndDebug(tree, 0);
+    addToRunTimeSelectionTable
+    (
+        evapotranspirationHeatTransferModel,
+        tree,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::scalar Foam::evapotranspirationHeatTransferModels::tree::Et() const
+{
+    return a_*q() + b_;
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::evapotranspirationHeatTransferModels::tree::leafArea
+(
+    const labelList& cells
+) const
+{
+    const volScalarField& LAD = getOrReadField(LADname_);
+
+    const scalarField& V = mesh().V();
+
+    auto tleafArea = tmp<scalarField>::New(cells.size(), Zero);
+    auto& leafArea = tleafArea.ref();
+
+    forAll(cells, i)
+    {
+        const label celli = cells[i];
+
+        leafArea[i] = LAD[celli]*V[celli];
+    }
+
+    return tleafArea;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::evapotranspirationHeatTransferModels::tree::tree
+(
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    evapotranspirationHeatTransferModel(dict, mesh),
+    a_(),
+    b_(),
+    lambda_(),
+    LADname_()
+{
+    Info<< "    Activating evapotranspiration heat transfer model: "
+        << typeName << endl;
+
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::scalarField>
+Foam::evapotranspirationHeatTransferModels::tree::Q
+(
+    const labelList& cells
+) const
+{
+    // Convert units from [MJ g/hr] to [J kg/s]
+    static const scalar unitConverter = scalar(1000)/scalar(3600);
+
+    return unitConverter*lambda_*Et()*leafArea(cells);
+}
+
+
+bool Foam::evapotranspirationHeatTransferModels::tree::read
+(
+    const dictionary& dict
+)
+{
+    if (!evapotranspirationHeatTransferModel::read(dict))
+    {
+        return false;
+    }
+
+    const auto range = scalarMinMax::ge(SMALL);
+
+    a_ = dict.getOrDefault<scalar>("a", 0.3622);
+    b_ = dict.getOrDefault<scalar>("b", 60.758);
+    lambda_ = dict.getCheckOrDefault<scalar>("lambda", 2.44, range);
+    LADname_ = dict.getOrDefault<word>("LAD", "LAD");
+
+    (void) getOrReadField(LADname_);
+
+    return true;
+}
+
+
+// ************************************************************************* //

src/atmosphericModels/fvOptions/evapotranspirationHeatTransfer/evapotranspirationHeatTransferModels/tree/tree.H

@@ -0,0 +1,219 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::evapotranspirationHeatTransferModels::tree
+
+Description
+    Applies sources on temperature (\c T - incompressible) or energy
+    (\c h/e - compressible) equation to incorporate evapotranspiration
+    heat-transfer effects from the specified tree canopy.
+
+    The heat-transfer is calculated based on the following relations.
+    The heat transfer is given by the first equation and the
+    evapotranspiration is calculated linearly proportional to absorbed
+    solar radiation (direct and diffusive) with the empirical relation
+    shown in the second equation below.
+
+    \f[
+        Q = \lambda E_t \frac{1000}{3600} \int_{cellZone} LAD \, dV
+    \f]
+
+    where
+    \vartable
+      Q         | Heat transfer from tree crown                 [W]
+      \lambda   | Specific latent heat of vaporisation          [MJ/kg]
+      E_t       | Area-averaged evapotranspiration              [g/m^2/hr]
+      LAD       | Leaf area density                             [m^2/m^3]
+      V         | Volume                                        [m^3]
+    \endvartable
+
+    with
+
+    \f[
+        E_t = a R_n + b
+    \f]
+
+    where
+    \vartable
+      R_n    | Heat flux through tree crown                    [W/m^2]
+      a      | Linear-regression coefficient - slope parameter [g/W/hr]
+      b      | Linear-regression coefficient - intercept parameter [g/m^2/hr]
+    \endvartable
+
+    Sources applied to either of the below, if exist:
+    \verbatim
+      T         | Temperature                                [K]
+      e         | Internal energy                            [m^2/s^2]
+      h         | Enthalphy                                  [m^2/s^2]
+    \endverbatim
+
+    Required fields:
+    \verbatim
+      T         | Temperature                                [K]
+      e         | Internal energy                            [m^2/s^2]
+      h         | Enthalphy                                  [m^2/s^2]
+      LAD       | Leaf area density                          [m^2/m^3]
+    \endverbatim
+
+Usage
+    Example by using \c constant/fvOptions:
+    \verbatim
+    evapotranspirationHeatTransfer1
+    {
+        // Inherited entries
+        ...
+
+        // Mandatory entries
+        model               tree;
+
+        // Optional entries
+        a                   <scalar>;
+        b                   <scalar>;
+        lambda              <scalar>;
+        LAD                 <word>;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property | Description                         | Type | Reqd | Deflt
+      model    | Model name: tree                    | word | yes  | -
+      a | Linear-regression coefficient - slope parameter [g/W/hr] | scalar <!--
+        -->                                                   | no | 0.3622
+      b | Linear-regression coefficient - intercept parameter [g/m^2/hr] <!--
+        -->                                         | scalar  | no | 60.758
+      lambda   | Specific latent heat of vaporisation [MJ/kg] | scalar <!--
+        -->                                                   | no | 2.44
+      LAD      | Name of leaf-area density field       | word | no | LAD
+    \endtable
+
+    The inherited entries are elaborated in:
+      - \link evapotranspirationHeatTransfer.H \endlink
+      - \link evapotranspirationHeatTransferModel.H \endlink
+
+SourceFiles
+    tree.C
+    treeTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_evapotranspirationHeatTransferModels_tree_H
+#define Foam_evapotranspirationHeatTransferModels_tree_H
+
+#include "evapotranspirationHeatTransferModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace evapotranspirationHeatTransferModels
+{
+
+/*---------------------------------------------------------------------------*\
+                            Class tree Declaration
+\*---------------------------------------------------------------------------*/
+
+class tree
+:
+    public evapotranspirationHeatTransferModel
+{
+    // Private Data
+
+        //- Linear-regression coefficient - slope parameter [g/W/hr]
+        scalar a_;
+
+        //- Linear-regression coefficient - intercept parameter [g/m^2/hr]
+        scalar b_;
+
+        //- Specific latent heat of vaporisation [MJ/kg]
+        scalar lambda_;
+
+        //- Name of leaf-area density field
+        word LADname_;
+
+
+    // Private Member Functions
+
+        //- Return area-averaged transpiration [g/m^2/hr]
+        //  Calculated from empirical regressions with heat flux
+        scalar Et() const;
+
+        //- Return volume weighted leaf area density (LAD) [m^2]
+        //  LAD [m^2/m^3]: one-sided area of leaves per unit volume
+        tmp<scalarField> leafArea(const labelList& cells) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("tree");
+
+
+    // Constructors
+
+        //- Construct from components
+        tree
+        (
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+        //- No copy construct
+        tree(const tree&) = delete;
+
+        //- No copy assignment
+        void operator=(const tree&) = delete;
+
+
+    //- Destructor
+    virtual ~tree() = default;
+
+
+    // Member Functions
+
+    // Evaluation
+
+        //- Return heat-transfer rate for speficied cells [J/s]
+        virtual tmp<scalarField> Q(const labelList& cells) const;
+
+
+    // I-O
+
+        //- Read the dictionary
+        virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace evapotranspirationHeatTransferModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/atmosphericModels/fvOptions/treeTurbulence/treeTurbulence.C

@@ -0,0 +1,234 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "treeTurbulence.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+    defineTypeNameAndDebug(treeTurbulence, 0);
+    addToRunTimeSelectionTable(option, treeTurbulence, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+Foam::volScalarField& Foam::fv::treeTurbulence::getOrReadField
+(
+    const word& fieldName
+) const
+{
+    auto* ptr = mesh_.getObjectPtr<volScalarField>(fieldName);
+
+    if (!ptr)
+    {
+        ptr = new volScalarField
+        (
+            IOobject
+            (
+                fieldName,
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::MUST_READ,
+                IOobject::AUTO_WRITE
+            ),
+            mesh_
+        );
+        mesh_.objectRegistry::store(ptr);
+    }
+
+    return *ptr;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fv::treeTurbulence::treeTurbulence
+(
+    const word& sourceName,
+    const word& modelType,
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    fv::cellSetOption(sourceName, modelType, dict, mesh),
+    isEpsilon_(true),
+    betaP_(),
+    betaD_(),
+    Ceps1_(),
+    Ceps2_(),
+    betaStar_(),
+    CdName_(),
+    LADname_()
+{
+    read(dict);
+
+    const auto* turbPtr = mesh_.findObject<turbulenceModel>
+    (
+        turbulenceModel::propertiesName
+    );
+
+    if (!turbPtr)
+    {
+        FatalErrorInFunction
+            << "Unable to find a turbulence model."
+            << abort(FatalError);
+    }
+
+    fieldNames_.resize(2);
+
+    tmp<volScalarField> tepsilon = turbPtr->epsilon();
+    tmp<volScalarField> tomega = turbPtr->omega();
+
+    if (!tepsilon.isTmp())
+    {
+        fieldNames_[0] = tepsilon().name();
+    }
+    else if (!tomega.isTmp())
+    {
+        isEpsilon_ = false;
+        fieldNames_[0] = tomega().name();
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Unable to find epsilon or omega field." << nl
+            << abort(FatalError);
+    }
+
+    fieldNames_[1] = turbPtr->k()().name();
+
+    fv::option::resetApplied();
+
+    Log << "    Applying treeTurbulence to: "
+        << fieldNames_[0] << " and " << fieldNames_[1]
+        << endl;
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::fv::treeTurbulence::addSup
+(
+    fvScalarMatrix& eqn,
+    const label fieldi
+)
+{
+    if (fieldi == 1)
+    {
+        kSource(geometricOneField(), geometricOneField(), eqn, fieldi);
+        return;
+    }
+
+    if (isEpsilon_)
+    {
+        epsilonSource(geometricOneField(), geometricOneField(), eqn, fieldi);
+    }
+    else
+    {
+        omegaSource(geometricOneField(), geometricOneField(), eqn, fieldi);
+    }
+}
+
+
+void Foam::fv::treeTurbulence::addSup
+(
+    const volScalarField& rho,
+    fvScalarMatrix& eqn,
+    const label fieldi
+)
+{
+    if (fieldi == 1)
+    {
+        kSource(geometricOneField(), rho, eqn, fieldi);
+        return;
+    }
+
+    if (isEpsilon_)
+    {
+        epsilonSource(geometricOneField(), rho, eqn, fieldi);
+    }
+    else
+    {
+        omegaSource(geometricOneField(), rho, eqn, fieldi);
+    }
+}
+
+
+void Foam::fv::treeTurbulence::addSup
+(
+    const volScalarField& alpha,
+    const volScalarField& rho,
+    fvScalarMatrix& eqn,
+    const label fieldi
+)
+{
+    if (fieldi == 1)
+    {
+        kSource(alpha, rho, eqn, fieldi);
+        return;
+    }
+
+    if (isEpsilon_)
+    {
+        epsilonSource(alpha, rho, eqn, fieldi);
+    }
+    else
+    {
+        omegaSource(alpha, rho, eqn, fieldi);
+    }
+}
+
+
+bool Foam::fv::treeTurbulence::read(const dictionary& dict)
+{
+    if (!fv::cellSetOption::read(dict))
+    {
+        return false;
+    }
+
+    betaP_ = dict.getOrDefault<scalar>("betaP", 1.0);
+    betaD_ = dict.getOrDefault<scalar>("betaD", 4.0);
+    Ceps1_ = dict.getOrDefault<scalar>("Ceps1", 0.9);
+    Ceps2_ = dict.getOrDefault<scalar>("Ceps2", 0.9);
+    betaStar_ = dict.getOrDefault<scalar>("betaStar", 0.09);
+    CdName_ = dict.getOrDefault<word>("Cd", "Cd");
+    LADname_ = dict.getOrDefault<word>("LAD", "LAD");
+
+    (void) getOrReadField(CdName_);
+    (void) getOrReadField(LADname_);
+
+    return true;
+}
+
+
+// ************************************************************************* //

src/atmosphericModels/fvOptions/treeTurbulence/treeTurbulence.H

@@ -0,0 +1,262 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::fv::treeTurbulence
+
+Group
+    grpFvOptionsSources
+
+Description
+    Applies sources on \c k and either \c epsilon or \c omega
+    to incorporate effects of trees on atmospheric boundary layer.
+
+    Corrections applied to:
+    \verbatim
+      k         | Turbulent kinetic energy                  [m^2/s^2]
+    \endverbatim
+
+    Corrections applied to either of the below, if exist:
+    \verbatim
+      epsilon   | Turbulent kinetic energy dissipation rate [m^2/s^3]
+      omega     | Specific dissipation rate                 [1/s]
+    \endverbatim
+
+    Required fields:
+    \verbatim
+      k         | Turbulent kinetic energy                  [m^2/s^2]
+      Cd        | Canopy drag coefficient                   [-]
+      LAD       | Leaf area density                         [m^2/m^3]
+      epsilon   | Turbulent kinetic energy dissipation rate [m^2/s^3]
+      omega     | Specific dissipation rate                 [1/s]
+    \endverbatim
+
+    References:
+    \verbatim
+        Governing equations (tag:BSG):
+            Brozovsky, J., Simonsen, A., & Gaitani, N. (2021).
+            Validation of a CFD model for the evaluation of urban microclimate
+            at high latitudes: A case study in Trondheim, Norway.
+            Building and Environment, 205, 108175.
+            DOI:10.1016/j.buildenv.2021.108175
+    \endverbatim
+
+Usage
+    Example by using \c constant/fvOptions:
+    \verbatim
+    treeTurbulence1
+    {
+        // Mandatory entries
+        type         treeTurbulence;
+
+        // Optional entries
+        Cd           <word>;
+        LAD          <word>;
+        betaP        <scalar>;
+        betaD        <scalar>;
+        Ceps1        <scalar>;
+        Ceps2        <scalar>;
+        betaStar     <scalar>;
+
+        // Inherited entries
+        ...
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property | Description                         | Type | Reqd | Deflt
+      type     | Type name: treeTurbulence           | word | yes  | -
+      Cd       | Name of operand canopy drag coefficient field | word | no | Cd
+      LAD      | Name of operand leaf area density field | word | no | LAD
+      betaP    | Model coefficient                 | scalar | no   | 1.0
+      betaD    | Model coefficient                 | scalar | no   | 4.0
+      Ceps1    | Model coefficient                 | scalar | no   | 0.9
+      Ceps2    | Model coefficient                 | scalar | no   | 0.9
+      betaStar | Model coefficient                 | scalar | no   | 0.09
+    \endtable
+
+    The inherited entries are elaborated in:
+      - \link fvOption.H \endlink
+      - \link cellSetOption.H \endlink
+
+SourceFiles
+    treeTurbulence.C
+    treeTurbulenceTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fv_treeTurbulence_H
+#define fv_treeTurbulence_H
+
+#include "cellSetOption.H"
+#include "turbulentTransportModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class treeTurbulence Declaration
+\*---------------------------------------------------------------------------*/
+
+class treeTurbulence
+:
+    public fv::cellSetOption
+{
+    // Private Data
+
+        //- Internal flag to determine the working field is epsilon or omega
+        bool isEpsilon_;
+
+        // Model coefficients
+
+            scalar betaP_;
+            scalar betaD_;
+            scalar Ceps1_;
+            scalar Ceps2_;
+            scalar betaStar_;
+
+        //- Name of operand canopy drag coefficient field
+        word CdName_;
+
+        //- Name of operand leaf area density field
+        word LADname_;
+
+
+     // Private Member Functions
+
+        //- Return requested field from the object registry
+        //- or read+register the field to the object registry
+        volScalarField& getOrReadField(const word& fieldName) const;
+
+        //- Apply sources to turbulent kinetic energy
+        template<class AlphaFieldType, class RhoFieldType>
+        void kSource
+        (
+            const AlphaFieldType& alpha,
+            const RhoFieldType& rho,
+            fvScalarMatrix& eqn,
+            const label fieldi
+        ) const;
+
+        //- Apply sources to turbulent kinetic energy dissipation rate
+        template<class AlphaFieldType, class RhoFieldType>
+        void epsilonSource
+        (
+            const AlphaFieldType& alpha,
+            const RhoFieldType& rho,
+            fvScalarMatrix& eqn,
+            const label fieldi
+        ) const;
+
+        //- Apply sources to specific dissipation rate
+        template<class AlphaFieldType, class RhoFieldType>
+        void omegaSource
+        (
+            const AlphaFieldType& alpha,
+            const RhoFieldType& rho,
+            fvScalarMatrix& eqn,
+            const label fieldi
+        ) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("treeTurbulence");
+
+
+    // Constructors
+
+        //- Construct from explicit source name and mesh
+        treeTurbulence
+        (
+            const word& sourceName,
+            const word& modelType,
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+        //- No copy construct
+        treeTurbulence(const treeTurbulence&) = delete;
+
+        //- No copy assignment
+        void operator=(const treeTurbulence&) = delete;
+
+
+    // Member Functions
+
+        //- Add explicit contribution to epsilon or omega equation
+        //- for incompressible flow computations
+         virtual void addSup
+        (
+            fvScalarMatrix& eqn,
+            const label fieldi
+        );
+
+        //- Add explicit contribution to epsilon or omega equation
+        //- for compressible flow computations
+        virtual void addSup
+        (
+            const volScalarField& rho,
+            fvScalarMatrix& eqn,
+            const label fieldi
+        );
+
+        //- Add explicit contribution to epsilon or omega equation
+        //- for multiphase flow computations
+        virtual void addSup
+        (
+            const volScalarField& alpha,
+            const volScalarField& rho,
+            fvScalarMatrix& eqn,
+            const label fieldi
+        );
+
+        //- Read source dictionary
+        virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "treeTurbulenceTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/atmosphericModels/fvOptions/treeTurbulence/treeTurbulenceTemplates.C

@@ -0,0 +1,124 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "treeTurbulence.H"
+#include "volFields.H"
+#include "fvmSup.H"
+
+// * * * * * * * * * * * * * * *  Member Functions * * * * * * * * * * * * * //
+
+template<class AlphaFieldType, class RhoFieldType>
+void Foam::fv::treeTurbulence::kSource
+(
+    const AlphaFieldType& alpha,
+    const RhoFieldType& rho,
+    fvScalarMatrix& eqn,
+    const label fieldi
+) const
+{
+    const auto* turbPtr = mesh_.findObject<turbulenceModel>
+    (
+        turbulenceModel::propertiesName
+    );
+
+    const volScalarField& k = turbPtr->k();
+    const volVectorField& U = turbPtr->U();
+    const volScalarField magU(mag(U));
+    const volScalarField& Cd = getOrReadField(CdName_);
+    const volScalarField& LAD = getOrReadField(LADname_);
+
+    // (BSG:Eq. 8)
+    eqn +=
+        alpha*rho*LAD*Cd*betaP_*pow3(magU)
+      - fvm::Sp(alpha*rho*LAD*Cd*betaD_*magU, k);
+}
+
+
+template<class AlphaFieldType, class RhoFieldType>
+void Foam::fv::treeTurbulence::epsilonSource
+(
+    const AlphaFieldType& alpha,
+    const RhoFieldType& rho,
+    fvScalarMatrix& eqn,
+    const label fieldi
+) const
+{
+    const auto* turbPtr = mesh_.findObject<turbulenceModel>
+    (
+        turbulenceModel::propertiesName
+    );
+
+    const volScalarField& epsilon = turbPtr->epsilon();
+    const volScalarField& k = turbPtr->k();
+    const volVectorField& U = turbPtr->U();
+    const volScalarField magU(mag(U));
+    const volScalarField& Cd = getOrReadField(CdName_);
+    const volScalarField& LAD = getOrReadField(LADname_);
+
+    // (BSG:Eq. 9)
+    eqn +=
+        fvm::Sp
+        (
+            alpha*rho*LAD*Cd*(Ceps1_*betaP_/k*pow3(magU) - Ceps2_*betaD_*magU),
+            epsilon
+        );
+}
+
+
+template<class AlphaFieldType, class RhoFieldType>
+void Foam::fv::treeTurbulence::omegaSource
+(
+    const AlphaFieldType& alpha,
+    const RhoFieldType& rho,
+    fvScalarMatrix& eqn,
+    const label fieldi
+) const
+{
+    const auto* turbPtr = mesh_.findObject<turbulenceModel>
+    (
+        turbulenceModel::propertiesName
+    );
+
+    const volScalarField& omega = turbPtr->omega();
+    const volScalarField& k = turbPtr->k();
+    const volVectorField& U = turbPtr->U();
+    const volScalarField magU(mag(U));
+    const volScalarField& Cd = getOrReadField(CdName_);
+    const volScalarField& LAD = getOrReadField(LADname_);
+
+    // (derived from BSG:Eq. 9 by assuming epsilon = betaStar*omega*k)
+    eqn +=
+        fvm::Sp
+        (
+            alpha*rho*LAD*Cd*betaStar_
+           *(Ceps1_*betaP_*pow3(magU) - Ceps2_*betaD_*k*magU),
+            omega
+        );
+}
+
+
+// ************************************************************************* //

src/dynamicMesh/Make/files

@@ -109,19 +109,6 @@ motionSolvers/displacement/solidBody/solidBodyMotionSolver.C
 motionSolvers/displacement/solidBody/multiSolidBodyMotionSolver.C
 motionSolvers/displacement/codedPoints0/codedPoints0MotionSolver.C
 
-solidBodyMotionFunctions = motionSolvers/displacement/solidBody/solidBodyMotionFunctions
-$(solidBodyMotionFunctions)/solidBodyMotionFunction/solidBodyMotionFunction.C
-$(solidBodyMotionFunctions)/solidBodyMotionFunction/solidBodyMotionFunctionNew.C
-$(solidBodyMotionFunctions)/SDA/SDA.C
-$(solidBodyMotionFunctions)/tabulated6DoFMotion/tabulated6DoFMotion.C
-$(solidBodyMotionFunctions)/linearMotion/linearMotion.C
-$(solidBodyMotionFunctions)/drivenLinearMotion/drivenLinearMotion.C
-$(solidBodyMotionFunctions)/rotatingMotion/rotatingMotion.C
-$(solidBodyMotionFunctions)/axisRotationMotion/axisRotationMotion.C
-$(solidBodyMotionFunctions)/multiMotion/multiMotion.C
-$(solidBodyMotionFunctions)/oscillatingLinearMotion/oscillatingLinearMotion.C
-$(solidBodyMotionFunctions)/oscillatingRotatingMotion/oscillatingRotatingMotion.C
-
 motionSolvers/displacement/solidBody/pointPatchFields/derived/solidBodyMotionDisplacement/solidBodyMotionDisplacementPointPatchVectorField.C
 
 createShellMesh/createShellMesh.C

src/dynamicMesh/polyMeshAdder/polyMeshAdder.C

@@ -2817,7 +2817,7 @@ void Foam::polyMeshAdder::add
                 if (newNei < newOwn)
                 {
                     std::swap(newOwn, newNei);
-                    newFace = newFace.reverseFace();
+                    newFace.flip();
                     flipFaceFlux = !flipFaceFlux;
                     for (bool& flip : flips)
                     {

src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.C

@@ -109,7 +109,21 @@ Foam::cyclicACMIFvPatchField<Type>::cyclicACMIFvPatchField
                 << exit(FatalIOError);
         }
 
+        // Tricky: avoid call to evaluate without call to initEvaluate.
+        // For now just disable the localConsistency to make it use the
+        // old logic (ultimately calls the fully self contained
+        // patchNeighbourField)
+
+        int& consistency =
+            GeometricField<Type, fvPatchField, volMesh>::
+            Boundary::localConsistency;
+
+        const int oldConsistency = consistency;
+        consistency = 0;
+
         this->evaluate(Pstream::commsTypes::blocking);
+
+        consistency = oldConsistency;
     }
 }
 

src/finiteVolume/fields/fvPatchFields/constraint/cyclicAMI/cyclicAMIFvPatchField.C

@@ -92,7 +92,21 @@ Foam::cyclicAMIFvPatchField<Type>::cyclicAMIFvPatchField
     {
         if (this->coupled())
         {
+            // Tricky: avoid call to evaluate without call to initEvaluate.
+            // For now just disable the localConsistency to make it use the
+            // old logic (ultimately calls the fully self contained
+            // patchNeighbourField)
+
+            int& consistency =
+                GeometricField<Type, fvPatchField, volMesh>::
+                Boundary::localConsistency;
+
+            const int oldConsistency = consistency;
+            consistency = 0;
+
             this->evaluate(UPstream::commsTypes::nonBlocking);
+
+            consistency = oldConsistency;
         }
         else
         {

src/meshTools/Make/files

@@ -354,4 +354,18 @@ tetOverlapVolume/tetOverlapVolume.C
 
 triangulatedPatch/triangulatedPatch.C
 
+solidBodyMotion = solidBodyMotionFunctions
+$(solidBodyMotion)/solidBodyMotionFunction/solidBodyMotionFunction.C
+$(solidBodyMotion)/solidBodyMotionFunction/solidBodyMotionFunctionNew.C
+$(solidBodyMotion)/SDA/SDA.C
+$(solidBodyMotion)/tabulated6DoFMotion/tabulated6DoFMotion.C
+$(solidBodyMotion)/linearMotion/linearMotion.C
+$(solidBodyMotion)/drivenLinearMotion/drivenLinearMotion.C
+$(solidBodyMotion)/rotatingMotion/rotatingMotion.C
+$(solidBodyMotion)/axisRotationMotion/axisRotationMotion.C
+$(solidBodyMotion)/multiMotion/multiMotion.C
+$(solidBodyMotion)/oscillatingLinearMotion/oscillatingLinearMotion.C
+$(solidBodyMotion)/oscillatingRotatingMotion/oscillatingRotatingMotion.C
+
+
 LIB = $(FOAM_LIBBIN)/libmeshTools

src/meshTools/Make/options

@@ -1,7 +1,6 @@
 EXE_INC = \
     -I$(LIB_SRC)/fileFormats/lnInclude \
-    -I$(LIB_SRC)/surfMesh/lnInclude \
-    -I$(LIB_SRC)/dynamicMesh/lnInclude
+    -I$(LIB_SRC)/surfMesh/lnInclude
 
 LIB_LIBS = \
     -lfileFormats \

src/meshTools/solidBodyMotionFunctions/solidBodyMotionFunction/solidBodyMotionFunction.H

@@ -39,12 +39,12 @@ Description
 
 SourceFiles
     solidBodyMotionFunction.C
-    dynamicFvMeshNew.C
+    solidBodyMotionFunctionNew.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef solidBodyMotionFunction_H
-#define solidBodyMotionFunction_H
+#ifndef Foam_solidBodyMotionFunction_H
+#define Foam_solidBodyMotionFunction_H
 
 #include "Time.H"
 #include "dictionary.H"
@@ -101,7 +101,7 @@ public:
 
     // Constructors
 
-        //- Construct from the SBMFCoeffs dictionary and Time
+        //- Construct from the coefficients dictionary and Time
         solidBodyMotionFunction
         (
             const dictionary& SBMFCoeffs,
@@ -114,10 +114,28 @@ public:
 
     // Selectors
 
-        //- Select constructed from the SBMFCoeffs dictionary and Time
+        //- Construct and dispatch motionType with dictionary and Time.
+        //  Returns nullptr if motionType is empty
         static autoPtr<solidBodyMotionFunction> New
         (
-            const dictionary& SBMFCoeffs,
+            const word& motionType,
+            const dictionary& dict,
+            const Time& runTime
+        );
+
+        //- Select "solidBodyMotionFunction" type from dictionary
+        //- and create with Time.
+        static autoPtr<solidBodyMotionFunction> New
+        (
+            const dictionary& dict,
+            const Time& runTime
+        );
+
+        //- Select "solidBodyMotionFunction" type (if present) from dictionary
+        //- and create with Time.
+        static autoPtr<solidBodyMotionFunction> NewIfPresent
+        (
+            const dictionary& dict,
             const Time& runTime
         );
 

src/meshTools/solidBodyMotionFunctions/solidBodyMotionFunction/solidBodyMotionFunctionNew.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,15 +30,18 @@ License
 
 // * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
 
-Foam::autoPtr<Foam::solidBodyMotionFunction> Foam::solidBodyMotionFunction::New
+Foam::autoPtr<Foam::solidBodyMotionFunction>
+Foam::solidBodyMotionFunction::New
 (
+    const word& motionType,
     const dictionary& dict,
     const Time& runTime
 )
 {
-    const word motionType(dict.get<word>("solidBodyMotionFunction"));
-
-    Info<< "Selecting solid-body motion function " << motionType << endl;
+    if (motionType.empty())
+    {
+        return nullptr;
+    }
 
     auto* ctorPtr = dictionaryConstructorTable(motionType);
 
@@ -53,8 +56,45 @@ Foam::autoPtr<Foam::solidBodyMotionFunction> Foam::solidBodyMotionFunction::New
         ) << exit(FatalIOError);
     }
 
+    Info<< "Selecting solid-body motion function " << motionType << endl;
+
     return autoPtr<solidBodyMotionFunction>(ctorPtr(dict, runTime));
 }
 
 
+Foam::autoPtr<Foam::solidBodyMotionFunction>
+Foam::solidBodyMotionFunction::New
+(
+    const dictionary& dict,
+    const Time& runTime
+)
+{
+    return New
+    (
+        dict.get<word>("solidBodyMotionFunction", keyType::LITERAL),
+        dict,
+        runTime
+    );
+}
+
+
+Foam::autoPtr<Foam::solidBodyMotionFunction>
+Foam::solidBodyMotionFunction::NewIfPresent
+(
+    const dictionary& dict,
+    const Time& runTime
+)
+{
+    word motionType;
+    dict.readIfPresent("solidBodyMotionFunction", motionType, keyType::LITERAL);
+
+    if (motionType.empty())
+    {
+        return nullptr;
+    }
+
+    return New(motionType, dict, runTime);
+}
+
+
 // ************************************************************************* //

src/meshTools/topoSet/cellSources/boundaryToCell/boundaryToCell.C

@@ -90,10 +90,10 @@ Foam::boundaryToCell::boundaryToCell(const polyMesh& mesh)
 Foam::boundaryToCell::boundaryToCell
 (
     const polyMesh& mesh,
-    const dictionary&
+    const dictionary& dict
 )
 :
-    topoSetCellSource(mesh)
+    topoSetCellSource(mesh, dict)
 {}
 
 

src/meshTools/topoSet/cellSources/boxToCell/boxToCell.C

@@ -92,7 +92,8 @@ static void readBoxDim(const dictionary& dict, treeBoundBox& bb)
 
 void Foam::boxToCell::combine(topoSet& set, const bool add) const
 {
-    const pointField& ctrs = mesh_.cellCentres();
+    const tmp<pointField> tctrs(this->transform(mesh_.cellCentres()));
+    const pointField& ctrs = tctrs();
 
     forAll(ctrs, elemi)
     {
@@ -138,7 +139,7 @@ Foam::boxToCell::boxToCell
     const dictionary& dict
 )
 :
-    topoSetCellSource(mesh),
+    topoSetCellSource(mesh, dict),
     bbs_()
 {
     // Accept 'boxes', 'box' or 'min/max'

src/meshTools/topoSet/cellSources/boxToCell/boxToCell.H

@@ -50,7 +50,7 @@ Usage
         source      boxToCell;
 
         // Conditional mandatory entries
-        // Select either of the below
+        // Select one of the below
 
         // Option-1
         boxes

src/meshTools/topoSet/cellSources/cellToCell/cellToCell.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -60,7 +60,8 @@ Foam::cellToCell::cellToCell
 )
 :
     topoSetCellSource(mesh),
-    names_(one{}, setName)
+    names_(Foam::one{}, setName),
+    isZone_(false)
 {}
 
 
@@ -70,16 +71,10 @@ Foam::cellToCell::cellToCell
     const dictionary& dict
 )
 :
-    topoSetCellSource(mesh),
-    names_()
-{
-    // Look for 'sets' or 'set'
-    if (!dict.readIfPresent("sets", names_))
-    {
-        names_.resize(1);
-        dict.readEntry("set", names_.front());
-    }
-}
+    topoSetCellSource(mesh, dict),
+    names_(),
+    isZone_(topoSetSource::readNames(dict, names_))
+{}
 
 
 Foam::cellToCell::cellToCell
@@ -89,7 +84,8 @@ Foam::cellToCell::cellToCell
 )
 :
     topoSetCellSource(mesh),
-    names_(one{}, word(checkIs(is)))
+    names_(Foam::one{}, word(checkIs(is))),
+    isZone_(false)
 {}
 
 
@@ -105,30 +101,44 @@ void Foam::cellToCell::applyToSet
     {
         if (verbose_)
         {
-            Info<< "    Adding all elements of cell sets: "
+            Info<< "    Adding all elements of cell "
+                << (isZone_ ? "zones: " : "sets: ")
                 << flatOutput(names_) << nl;
         }
 
         for (const word& setName : names_)
         {
-            cellSet loadedSet(mesh_, setName);
-
-            set.addSet(loadedSet);
+            if (isZone_)
+            {
+                set.addSet(mesh_.cellZones()[setName]);
+            }
+            else
+            {
+                cellSet loadedSet(mesh_, setName);
+                set.addSet(loadedSet);
+            }
         }
     }
     else if (action == topoSetSource::SUBTRACT)
     {
         if (verbose_)
         {
-            Info<< "    Removing all elements of cell sets: "
+            Info<< "    Removing all elements of cell "
+                << (isZone_ ? "zones: " : "sets: ")
                 << flatOutput(names_) << nl;
         }
 
         for (const word& setName : names_)
         {
-            cellSet loadedSet(mesh_, setName);
-
-            set.subtractSet(loadedSet);
+            if (isZone_)
+            {
+                set.subtractSet(mesh_.cellZones()[setName]);
+            }
+            else
+            {
+                cellSet loadedSet(mesh_, setName);
+                set.subtractSet(loadedSet);
+            }
         }
     }
 }

src/meshTools/topoSet/cellSources/cellToCell/cellToCell.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,7 +28,8 @@ Class
     Foam::cellToCell
 
 Description
-    A \c topoSetCellSource to select all the cells from given \c cellSet(s).
+    A \c topoSetCellSource to select all the cells from given \c cellSet(s)
+    or \c cellZone(s).
 
     Operands:
     \table
@@ -50,7 +51,7 @@ Usage
         source      cellToCell;
 
         // Conditional mandatory entries
-        // Select either of the below
+        // Select one of the below
 
         // Option-1
         sets
@@ -61,7 +62,18 @@ Usage
         );
 
         // Option-2
-        set <cellSetName>;
+        zones
+        (
+            <cellZoneName0>
+            <cellZoneName1>
+            ...
+        );
+
+        // Option-3
+        set     <cellSetName>;
+
+        // Option-4
+        zone    <cellZoneName>;
     }
     \endverbatim
 
@@ -81,17 +93,15 @@ Usage
       subtract | Remove selected cells from this cellSet
     \endverbatim
 
-    Options for the conditional mandatory entries:
+    Options for the conditional mandatory entries (in order of precedence):
     \verbatim
       Entry    | Description                    | Type     | Req'd  | Dflt
       sets     | Names of input cellSets        | wordList | cond'l | -
+      zones    | Names of input cellZones       | wordList | cond'l | -
       set      | Name of input cellSet          | word     | cond'l | -
+      zone     | Name of input cellZone         | word     | cond'l | -
     \endverbatim
 
-Note
-    The order of precedence among the conditional mandatory entries from the
-    highest to the lowest is \c sets, and \c set.
-
 See also
     - Foam::topoSetSource
     - Foam::topoSetCellSource
@@ -101,8 +111,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef cellToCell_H
-#define cellToCell_H
+#ifndef Foam_cellToCell_H
+#define Foam_cellToCell_H
 
 #include "topoSetCellSource.H"
 
@@ -124,9 +134,12 @@ class cellToCell
         //- Add usage string
         static addToUsageTable usage_;
 
-        //- Names of sets to use
+        //- Names of sets or zones to use
         wordList names_;
 
+        //- Is name a set or a zone
+        const bool isZone_;
+
 
 public:
 

src/meshTools/topoSet/cellSources/clipPlaneToCell/clipPlaneToCell.C

@@ -69,7 +69,8 @@ void Foam::clipPlaneToCell::combine(topoSet& set, const bool add) const
 {
     // Cell centres above the plane
 
-    const pointField& ctrs = mesh_.cellCentres();
+    const tmp<pointField> tctrs(this->transform(mesh_.cellCentres()));
+    const pointField& ctrs = tctrs();
 
     forAll(ctrs, elemi)
     {
@@ -102,12 +103,9 @@ Foam::clipPlaneToCell::clipPlaneToCell
     const dictionary& dict
 )
 :
-    clipPlaneToCell
-    (
-        mesh,
-        dict.get<vector>("point"),
-        dict.get<vector>("normal")
-    )
+    topoSetCellSource(mesh, dict),
+    point_(dict.get<vector>("point")),
+    normal_(dict.get<vector>("normal"))
 {}
 
 

src/meshTools/topoSet/cellSources/cylinderToCell/cylinderToCell.C

@@ -68,7 +68,8 @@ Foam::topoSetSource::addToUsageTable Foam::cylinderToCell::usage_
 
 void Foam::cylinderToCell::combine(topoSet& set, const bool add) const
 {
-    const pointField& ctrs = mesh_.cellCentres();
+    const tmp<pointField> tctrs(this->transform(mesh_.cellCentres()));
+    const pointField& ctrs = tctrs();
 
     const vector axis = (point2_ - point1_);
     const scalar magAxis2 = magSqr(axis);
@@ -127,12 +128,12 @@ Foam::cylinderToCell::cylinderToCell
     const dictionary& dict
 )
 :
-    cylinderToCell
+    topoSetCellSource(mesh, dict),
+    point1_(dict.getCompat<point>("point1", {{"p1", -2112}})),
+    point2_(dict.getCompat<point>("point2", {{"p2", -2112}})),
+    radius_(dict.getCompat<scalar>("radius", {{"outerRadius", -2112}})),
+    innerRadius_
     (
-        mesh,
-        dict.getCompat<point>("point1", {{"p1", -2112}}),
-        dict.getCompat<point>("point2", {{"p2", -2112}}),
-        dict.getCompat<scalar>("radius", {{"outerRadius", -2112}}),
         dict.getCheckOrDefault<scalar>("innerRadius", 0, scalarMinMax::ge(0))
     )
 {}

src/meshTools/topoSet/cellSources/faceToCell/faceToCell.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -66,18 +66,14 @@ Foam::faceToCell::faceActionNames_
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::faceToCell::combine
+template<class Selector>
+void Foam::faceToCell::combineImpl
 (
     topoSet& set,
     const bool add,
-    const word& setName
+    const Selector& faceLabels
 ) const
 {
-    // Load the set
-    faceSet loadedSet(mesh_, setName);
-
-    const labelHashSet& faceLabels = loadedSet;
-
     // Handle owner/neighbour/any selection
     for (const label facei : faceLabels)
     {
@@ -104,7 +100,7 @@ void Foam::faceToCell::combine
     {
         // Count number of selected faces per cell.
 
-        Map<label> facesPerCell(loadedSet.size());
+        Map<label> facesPerCell(faceLabels.size());
 
         for (const label facei : faceLabels)
         {
@@ -134,6 +130,31 @@ void Foam::faceToCell::combine
 }
 
 
+void Foam::faceToCell::combine
+(
+    topoSet& set,
+    const bool add,
+    const word& setName
+) const
+{
+    if (isZone_)
+    {
+        const labelList& faceLabels = mesh_.faceZones()[setName];
+
+        combineImpl(set, add, faceLabels);
+    }
+    else
+    {
+        // Load the set
+        faceSet loadedSet(mesh_, setName);
+
+        const labelHashSet& faceLabels = loadedSet;
+
+        combineImpl(set, add, faceLabels);
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::faceToCell::faceToCell
@@ -144,7 +165,8 @@ Foam::faceToCell::faceToCell
 )
 :
     topoSetCellSource(mesh),
-    names_(one{}, setName),
+    names_(Foam::one{}, setName),
+    isZone_(false),
     option_(option)
 {}
 
@@ -155,17 +177,11 @@ Foam::faceToCell::faceToCell
     const dictionary& dict
 )
 :
-    topoSetCellSource(mesh),
+    topoSetCellSource(mesh, dict),
     names_(),
+    isZone_(topoSetSource::readNames(dict, names_)),
     option_(faceActionNames_.get("option", dict))
-{
-    // Look for 'sets' or 'set'
-    if (!dict.readIfPresent("sets", names_))
-    {
-        names_.resize(1);
-        dict.readEntry("set", names_.front());
-    }
-}
+{}
 
 
 Foam::faceToCell::faceToCell
@@ -175,7 +191,8 @@ Foam::faceToCell::faceToCell
 )
 :
     topoSetCellSource(mesh),
-    names_(one{}, word(checkIs(is))),
+    names_(Foam::one{}, word(checkIs(is))),
+    isZone_(false),
     option_(faceActionNames_.read(checkIs(is)))
 {}
 
@@ -192,8 +209,9 @@ void Foam::faceToCell::applyToSet
     {
         if (verbose_)
         {
-            Info<< "    Adding cells according to face sets: "
-                << flatOutput(names_) << endl;
+            Info<< "    Adding cells according to face "
+                << (isZone_ ? "zones: " : "sets: ")
+                << flatOutput(names_) << nl;
         }
 
         for (const word& setName : names_)
@@ -205,8 +223,9 @@ void Foam::faceToCell::applyToSet
     {
         if (verbose_)
         {
-            Info<< "    Removing cells according to face sets: "
-                << flatOutput(names_) << endl;
+            Info<< "    Removing cells according to face "
+                << (isZone_ ? "zones: " : "sets: ")
+                << flatOutput(names_) << nl;
         }
 
         for (const word& setName : names_)

src/meshTools/topoSet/cellSources/faceToCell/faceToCell.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -53,7 +53,7 @@ Usage
         option      <option>;
 
         // Conditional mandatory entries
-        // Select either of the below
+        // Select one of the below
 
         // Option-1
         sets
@@ -64,7 +64,18 @@ Usage
         );
 
         // Option-2
+        zones
+        (
+            <faceZoneName0>
+            <faceZoneName1>
+            ...
+        );
+
+        // Option-3
         set     <faceSetName>;
+
+        // Option-4
+        zone    <faceZoneName>;
     }
     \endverbatim
 
@@ -93,17 +104,15 @@ Usage
       neighbour | Cells that are neighbour of given faces
     \endverbatim
 
-    Options for the conditional mandatory entries:
+    Options for the conditional mandatory entries (in order of precedence):
     \verbatim
       Entry    | Description                | Type     | Req'd  | Dflt
       sets     | Names of input faceSets    | wordList | cond'l | -
+      zones    | Names of input faceZones   | wordList | cond'l | -
       set      | Name of input faceSet      | word     | cond'l | -
+      zone     | Name of input faceZone     | word     | cond'l | -
     \endverbatim
 
-Note
-    The order of precedence among the conditional mandatory entries from the
-    highest to the lowest is \c sets, and \c set.
-
 See also
     - Foam::topoSetSource
     - Foam::topoSetCellSource
@@ -114,8 +123,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef faceToCell_H
-#define faceToCell_H
+#ifndef Foam_faceToCell_H
+#define Foam_faceToCell_H
 
 #include "topoSetCellSource.H"
 #include "Enum.H"
@@ -153,15 +162,27 @@ private:
         //- Add usage string
         static addToUsageTable usage_;
 
-        //- Names of sets to use
+        //- Names of sets or zones to use
         wordList names_;
 
+        //- Is name a set or a zone
+        const bool isZone_;
+
         //- Option
         faceAction option_;
 
 
     // Private Member Functions
 
+        //- Depending on face to cell option add to or delete from cellSet.
+        template<class Selector>
+        void combineImpl
+        (
+            topoSet& set,
+            const bool add,
+            const Selector& faceLabels
+        ) const;
+
         //- Depending on face to cell option add to or delete from cellSet.
         void combine(topoSet& set, const bool add, const word& setName) const;
 

src/meshTools/topoSet/cellSources/faceZoneToCell/faceZoneToCell.C

@@ -214,7 +214,7 @@ Foam::faceZoneToCell::faceZoneToCell
     const dictionary& dict
 )
 :
-    topoSetCellSource(mesh),
+    topoSetCellSource(mesh, dict),
     zoneMatcher_(),
     option_(faceActionNames_.get("option", dict))
 {

src/meshTools/topoSet/cellSources/faceZoneToCell/faceZoneToCell.H

@@ -52,7 +52,7 @@ Usage
         option      <option>;
 
         // Conditional mandatory entries
-        // Select either of the below
+        // Select one of the below
 
         // Option-1
         zones
@@ -63,7 +63,7 @@ Usage
         );
 
         // Option-2
-        set <faceZoneName>;
+        zone <faceZoneName>;
     }
     \endverbatim
 

src/meshTools/topoSet/cellSources/fieldToCell/fieldToCell.C

@@ -148,13 +148,10 @@ Foam::fieldToCell::fieldToCell
     const dictionary& dict
 )
 :
-    fieldToCell
-    (
-        mesh,
-        dict.get<word>("field"),
-        dict.get<scalar>("min"),
-        dict.get<scalar>("max")
-    )
+    topoSetCellSource(mesh, dict),
+    fieldName_(dict.get<word>("field")),
+    min_(dict.get<scalar>("min")),
+    max_(dict.get<scalar>("max"))
 {}
 
 

src/meshTools/topoSet/cellSources/haloToCell/haloToCell.C

@@ -198,7 +198,8 @@ Foam::haloToCell::haloToCell
     const dictionary& dict
 )
 :
-    haloToCell(mesh, dict.getOrDefault<label>("steps", 1))
+    topoSetCellSource(mesh, dict),
+    steps_(dict.getOrDefault<label>("steps", 1))
 {}
 
 

src/meshTools/topoSet/cellSources/labelToCell/labelToCell.C

@@ -94,7 +94,8 @@ Foam::labelToCell::labelToCell
     const dictionary& dict
 )
 :
-    labelToCell(mesh, dict.get<labelList>("value"))
+    topoSetCellSource(mesh, dict),
+    labels_(dict.get<labelList>("value"))
 {}
 
 

src/meshTools/topoSet/cellSources/nbrToCell/nbrToCell.C

@@ -74,6 +74,7 @@ void Foam::nbrToCell::combine(topoSet& set, const bool add) const
     }
 
     const cellList& cells = mesh().cells();
+
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
     boolList isCoupled(mesh_.nBoundaryFaces(), false);
@@ -136,7 +137,8 @@ Foam::nbrToCell::nbrToCell
     const dictionary& dict
 )
 :
-    nbrToCell(mesh, dict.getCheck<label>("neighbours", labelMinMax::ge(1)))
+    topoSetCellSource(mesh, dict),
+    minNbrs_(dict.getCheck<label>("neighbours", labelMinMax::ge(1)))
 {}
 
 

src/meshTools/topoSet/cellSources/nearestToCell/nearestToCell.C

@@ -127,11 +127,8 @@ Foam::nearestToCell::nearestToCell
     const dictionary& dict
 )
 :
-    nearestToCell
-    (
-        mesh,
-        dict.get<pointField>("points")
-    )
+    topoSetCellSource(mesh, dict),
+    points_(dict.get<pointField>("points"))
 {}
 
 

src/meshTools/topoSet/cellSources/patchToCell/patchToCell.C

@@ -125,7 +125,7 @@ Foam::patchToCell::patchToCell
     const dictionary& dict
 )
 :
-    topoSetCellSource(mesh),
+    topoSetCellSource(mesh, dict),
     selectedPatches_()
 {
     // Look for 'patches' and 'patch', but accept 'name' as well

src/meshTools/topoSet/cellSources/patchToCell/patchToCell.H

@@ -48,7 +48,7 @@ Usage
         source      patchToCell;
 
         // Conditional mandatory entries
-        // Select either of the below
+        // Select one of the below
 
         // Option-1
         patches

src/meshTools/topoSet/cellSources/pointToCell/pointToCell.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -64,18 +64,14 @@ Foam::pointToCell::pointActionNames_
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::pointToCell::combine
+template<class Selector>
+void Foam::pointToCell::combineImpl
 (
     topoSet& set,
     const bool add,
-    const word& setName
+    const Selector& pointLabels
 ) const
 {
-    // Load the set
-    pointSet loadedSet(mesh_, setName);
-
-    const labelHashSet& pointLabels = loadedSet;
-
     // Handle any selection
     if (option_ == ANY)
     {
@@ -114,6 +110,31 @@ void Foam::pointToCell::combine
 }
 
 
+void Foam::pointToCell::combine
+(
+    topoSet& set,
+    const bool add,
+    const word& setName
+) const
+{
+    if (isZone_)
+    {
+        const labelList& pointLabels = mesh_.pointZones()[setName];
+
+        combineImpl(set, add, pointLabels);
+    }
+    else
+    {
+        // Load the set
+        pointSet loadedSet(mesh_, setName);
+
+        const labelHashSet& pointLabels = loadedSet;
+
+        combineImpl(set, add, pointLabels);
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::pointToCell::pointToCell
@@ -124,7 +145,8 @@ Foam::pointToCell::pointToCell
 )
 :
     topoSetCellSource(mesh),
-    names_(one{}, setName),
+    names_(Foam::one{}, setName),
+    isZone_(false),
     option_(option)
 {}
 
@@ -135,17 +157,11 @@ Foam::pointToCell::pointToCell
     const dictionary& dict
 )
 :
-    topoSetCellSource(mesh),
+    topoSetCellSource(mesh, dict),
     names_(),
+    isZone_(topoSetSource::readNames(dict, names_)),
     option_(pointActionNames_.get("option", dict))
-{
-    // Look for 'sets' or 'set'
-    if (!dict.readIfPresent("sets", names_))
-    {
-        names_.resize(1);
-        dict.readEntry("set", names_.front());
-    }
-}
+{}
 
 
 Foam::pointToCell::pointToCell
@@ -155,7 +171,8 @@ Foam::pointToCell::pointToCell
 )
 :
     topoSetCellSource(mesh),
-    names_(one{}, word(checkIs(is))),
+    names_(Foam::one{}, word(checkIs(is))),
+    isZone_(false),
     option_(pointActionNames_.read(checkIs(is)))
 {}
 
@@ -172,7 +189,8 @@ void Foam::pointToCell::applyToSet
     {
         if (verbose_)
         {
-            Info<< "    Adding cells according to point sets: "
+            Info<< "    Adding cells according to point "
+                << (isZone_ ? "zones: " : "sets: ")
                 << flatOutput(names_) << nl;
         }
 
@@ -185,7 +203,8 @@ void Foam::pointToCell::applyToSet
     {
         if (verbose_)
         {
-            Info<< "    Removing cells according to point sets: "
+            Info<< "    Removing cells according to point "
+                << (isZone_ ? "zones: " : "sets: ")
                 << flatOutput(names_) << nl;
         }
 

src/meshTools/topoSet/cellSources/pointToCell/pointToCell.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2012 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -52,7 +52,7 @@ Usage
         option      <option>;
 
         // Conditional mandatory entries
-        // Select either of the below
+        // Select one of the below
 
         // Option-1
         sets
@@ -62,8 +62,19 @@ Usage
             ...
         );
 
-        // Option-2
-        set   <pointSetName>;
+        // Option-3
+        zones
+        (
+            <pointZoneName0>
+            <pointZoneName1>
+            ...
+        );
+
+        // Option-3
+        set     <pointSetName>;
+
+        // Option-4
+        zone    <pointZoneName>;
     }
     \endverbatim
 
@@ -90,17 +101,15 @@ Usage
       edge    | Cells using an edge with both points in pointSet
     \endverbatim
 
-    Options for the conditional mandatory entries:
+    Options for the conditional mandatory entries (in order of precedence):
     \verbatim
       Entry    | Description                | Type     | Req'd  | Dflt
       sets     | Names of input pointSets   | wordList | cond'l | -
+      zones    | Names of input pointZones  | wordList | cond'l | -
       set      | Name of input pointSet     | word     | cond'l | -
+      zone     | Name of input pointZone    | word     | cond'l | -
     \endverbatim
 
-Note
-    The order of precedence among the conditional mandatory entries from the
-    highest to the lowest is \c sets, and \c set.
-
 See also
     - Foam::topoSetSource
     - Foam::topoSetCellSource
@@ -110,8 +119,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef pointToCell_H
-#define pointToCell_H
+#ifndef Foam_pointToCell_H
+#define Foam_pointToCell_H
 
 #include "topoSetCellSource.H"
 #include "Enum.H"
@@ -148,15 +157,27 @@ private:
 
         static const Enum<pointAction> pointActionNames_;
 
-        //- Names of sets to use
+        //- Names of sets or zones to use
         wordList names_;
 
+        //- Is name a set or a zone
+        const bool isZone_;
+
         //- Selection type
         pointAction option_;
 
 
     // Private Member Functions
 
+        //- Depending on point to cell option add to or delete from cellSet.
+        template<class Selector>
+        void combineImpl
+        (
+            topoSet& set,
+            const bool add,
+            const Selector& pointLabels
+        ) const;
+
         //- Depending on point-to-cell option add to or delete from cellSet.
         void combine(topoSet& set, const bool add, const word& setName) const;
 

src/meshTools/topoSet/cellSources/regionToCell/regionToCell.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2020 OpenCFD Ltd.
+    Copyright (C) 2015-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -346,18 +346,33 @@ void Foam::regionToCell::combine(topoSet& set, const bool add) const
     // Note: wip. Select cells first
     boolList selectedCell(mesh_.nCells(), true);
 
-    if (setName_.size() && setName_ != "none")
+    if (!setName_.empty() && setName_ != "none")
     {
-        Info<< "    Loading subset " << setName_
-            << " to delimit search region."
-            << endl;
-
-        cellSet subSet(mesh_, setName_);
-
-        selectedCell = false;
-        for (const label celli : subSet)
+        if (isZone_)
         {
-            selectedCell[celli] = true;
+            Info<< "    Using cellZone " << setName_
+                << " to delimit search region."
+                << endl;
+
+            selectedCell = false;
+            for (const label celli : mesh_.cellZones()[setName_])
+            {
+                selectedCell[celli] = true;
+            }
+        }
+        else
+        {
+            Info<< "    Loading cellSet " << setName_
+                << " to delimit search region."
+                << endl;
+
+            cellSet subSet(mesh_, setName_);
+
+            selectedCell = false;
+            for (const label celli : subSet)
+            {
+                selectedCell[celli] = true;
+            }
         }
     }
 
@@ -392,6 +407,7 @@ Foam::regionToCell::regionToCell
 :
     topoSetCellSource(mesh),
     setName_(setName),
+    isZone_(false),
     insidePoints_(insidePoints),
     nErode_(nErode)
 {}
@@ -403,14 +419,25 @@ Foam::regionToCell::regionToCell
     const dictionary& dict
 )
 :
-    topoSetCellSource(mesh),
-    setName_(dict.getOrDefault<word>("set", "none")),
+    topoSetCellSource(mesh, dict),
+    setName_(),
+    isZone_(false),
     insidePoints_
     (
         dict.getCompat<pointField>("insidePoints", {{ "insidePoint", 0 }})
     ),
     nErode_(dict.getCheckOrDefault<label>("nErode", 0, labelMinMax::ge(0)))
-{}
+{
+    // A single set or zone only!
+    if (dict.readIfPresent("set", setName_))
+    {
+        isZone_ = false;
+    }
+    else if (dict.readIfPresent("zone", setName_))
+    {
+        isZone_ = true;
+    }
+}
 
 
 Foam::regionToCell::regionToCell
@@ -421,6 +448,7 @@ Foam::regionToCell::regionToCell
 :
     topoSetCellSource(mesh),
     setName_(checkIs(is)),
+    isZone_(false),
     insidePoints_(checkIs(is)),
     nErode_(readLabel(checkIs(is)))
 {}

src/meshTools/topoSet/cellSources/regionToCell/regionToCell.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2012 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -59,6 +59,7 @@ Usage
 
         // Optional entries
         set       <cellSetName>;
+        zone      <cellZoneName>;
         nErode    <label>;
     }
     \endverbatim
@@ -93,8 +94,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef regionToCell_H
-#define regionToCell_H
+#ifndef Foam_regionToCell_H
+#define Foam_regionToCell_H
 
 #include "topoSetCellSource.H"
 #include "boolList.H"
@@ -119,8 +120,11 @@ class regionToCell
         //- Add usage string
         static addToUsageTable usage_;
 
-        //- Name of cellSet to keep to
-        const word setName_;
+        //- Name of cellSet or cellZone to keep to
+        word setName_;
+
+        //- Is name a set or a zone
+        bool isZone_;
 
         //- Coordinate(s) that is inside connected region
         const pointField insidePoints_;

src/meshTools/topoSet/cellSources/rotatedBoxToCell/rotatedBoxToCell.C

@@ -96,7 +96,8 @@ void Foam::rotatedBoxToCell::combine(topoSet& set, const bool add) const
 
     // Check whether cell centre is inside all faces of box.
 
-    const pointField& ctrs = mesh_.cellCentres();
+    const tmp<pointField> tctrs(this->transform(mesh_.cellCentres()));
+    const pointField& ctrs = tctrs();
 
     forAll(ctrs, celli)
     {
@@ -146,14 +147,11 @@ Foam::rotatedBoxToCell::rotatedBoxToCell
     const dictionary& dict
 )
 :
-    rotatedBoxToCell
-    (
-        mesh,
-        dict.get<point>("origin"),
-        dict.get<vector>("i"),
-        dict.get<vector>("j"),
-        dict.get<vector>("k")
-    )
+    topoSetCellSource(mesh, dict),
+    origin_(dict.get<point>("origin")),
+    i_(dict.get<vector>("i")),
+    j_(dict.get<vector>("j")),
+    k_(dict.get<vector>("k"))
 {}
 
 

src/meshTools/topoSet/cellSources/searchableSurfaceToCell/searchableSurfaceToCell.C

@@ -102,9 +102,13 @@ void Foam::searchableSurfaceToCell::combine(topoSet& set, const bool add) const
     {
         return;
     }
-    const pointField& ctrs = mesh_.cellCentres();
+
+    const tmp<pointField> tctrs(this->transform(mesh_.cellCentres()));
+    const pointField& ctrs = tctrs();
+
     const searchableSurface& s = *surf_;
 
+
     // Cell centres within the enclosing volumes
 
     List<volumeType> volTypes;
@@ -130,7 +134,7 @@ Foam::searchableSurfaceToCell::searchableSurfaceToCell
     const dictionary& dict
 )
 :
-    topoSetCellSource(mesh),
+    topoSetCellSource(mesh, dict),
     surf_
     (
         searchableSurface::New

src/meshTools/topoSet/cellSources/shapeToCell/shapeToCell.C

@@ -116,7 +116,8 @@ Foam::shapeToCell::shapeToCell
     const dictionary& dict
 )
 :
-    shapeToCell(mesh, dict.getCompat<word>("shape", {{"type", 1806}}))
+    topoSetCellSource(mesh, dict),
+    shape_(dict.getCompat<word>("shape", {{"type", 1806}}))
 {}
 
 

src/meshTools/topoSet/cellSources/sphereToCell/sphereToCell.C

@@ -68,7 +68,8 @@ Foam::topoSetSource::addToUsageTable Foam::sphereToCell::usage_
 
 void Foam::sphereToCell::combine(topoSet& set, const bool add) const
 {
-    const pointField& ctrs = mesh_.cellCentres();
+    const tmp<pointField> tctrs(this->transform(mesh_.cellCentres()));
+    const pointField& ctrs = tctrs();
 
     const scalar orad2 = sqr(radius_);
     const scalar irad2 = innerRadius_ > 0 ? sqr(innerRadius_) : -1;
@@ -110,11 +111,11 @@ Foam::sphereToCell::sphereToCell
     const dictionary& dict
 )
 :
-    sphereToCell
+    topoSetCellSource(mesh, dict),
+    origin_(dict.getCompat<vector>("origin", {{"centre", -1806}})),
+    radius_(dict.getCheck<scalar>("radius", scalarMinMax::ge(0))),
+    innerRadius_
     (
-        mesh,
-        dict.getCompat<vector>("origin", {{"centre", -1806}}),
-        dict.getCheck<scalar>("radius", scalarMinMax::ge(0)),
         dict.getCheckOrDefault<scalar>("innerRadius", 0, scalarMinMax::ge(0))
     )
 {}

src/meshTools/topoSet/cellSources/surfaceToCell/surfaceToCell.C

@@ -434,7 +434,7 @@ Foam::surfaceToCell::surfaceToCell
     const dictionary& dict
 )
 :
-    topoSetCellSource(mesh),
+    topoSetCellSource(mesh, dict),
     surfName_(dict.get<fileName>("file").expand()),
     outsidePoints_(dict.get<pointField>("outsidePoints")),
     includeCut_(dict.get<bool>("includeCut")),

src/meshTools/topoSet/cellSources/targetVolumeToCell/targetVolumeToCell.C

@@ -301,13 +301,10 @@ Foam::targetVolumeToCell::targetVolumeToCell
     const dictionary& dict
 )
 :
-    targetVolumeToCell
-    (
-        mesh,
-        dict.getCheck<scalar>("volume", scalarMinMax::ge(0)),
-        dict.get<vector>("normal"),
-        dict.getOrDefault<word>("set", "")
-    )
+    topoSetCellSource(mesh, dict),
+    vol_(dict.getCheck<scalar>("volume", scalarMinMax::ge(0))),
+    normal_(dict.get<vector>("normal")),
+    maskSetName_(dict.getOrDefault<word>("set", ""))
 {}
 
 

src/meshTools/topoSet/cellSources/zoneToCell/zoneToCell.C

@@ -181,7 +181,7 @@ Foam::zoneToCell::zoneToCell
     const dictionary& dict
 )
 :
-    topoSetCellSource(mesh),
+    topoSetCellSource(mesh, dict),
     zoneMatcher_()
 {
     // Look for 'zones' and 'zone', but accept 'name' as well

src/meshTools/topoSet/cellSources/zoneToCell/zoneToCell.H

@@ -50,7 +50,7 @@ Usage
         source      zoneToCell;
 
         // Conditional mandatory entries
-        // Select either of the below
+        // Select one of the below
 
         // Option-1
         zones

src/meshTools/topoSet/cellZoneSources/setToCellZone/setToCellZone.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -71,7 +71,7 @@ Foam::setToCellZone::setToCellZone
     const dictionary& dict
 )
 :
-    topoSetCellZoneSource(mesh),
+    topoSetCellZoneSource(mesh, dict),
     setName_(dict.get<word>("set"))
 {}
 
@@ -123,7 +123,7 @@ void Foam::setToCellZone::applyToSet
             {
                 if (!zoneSet.found(celli))
                 {
-                    newAddressing.append(celli);
+                    newAddressing.push_back(celli);
                 }
             }
 
@@ -144,11 +144,11 @@ void Foam::setToCellZone::applyToSet
             // Start off empty
             DynamicList<label> newAddressing(zoneSet.addressing().size());
 
-            forAll(zoneSet.addressing(), i)
+            for (const label celli : zoneSet.addressing())
             {
-                if (!loadedSet.found(zoneSet.addressing()[i]))
+                if (!loadedSet.found(celli))
                 {
-                    newAddressing.append(zoneSet.addressing()[i]);
+                    newAddressing.push_back(celli);
                 }
             }
             zoneSet.addressing().transfer(newAddressing);

src/meshTools/topoSet/faceSources/boundaryToFace/boundaryToFace.C

@@ -91,10 +91,10 @@ Foam::boundaryToFace::boundaryToFace(const polyMesh& mesh)
 Foam::boundaryToFace::boundaryToFace
 (
     const polyMesh& mesh,
-    const dictionary&
+    const dictionary& dict
 )
 :
-    topoSetFaceSource(mesh)
+    topoSetFaceSource(mesh, dict)
 {}
 
 

src/meshTools/topoSet/faceSources/boxToFace/boxToFace.C

@@ -92,7 +92,8 @@ static void readBoxDim(const dictionary& dict, treeBoundBox& bb)
 
 void Foam::boxToFace::combine(topoSet& set, const bool add) const
 {
-    const pointField& ctrs = mesh_.faceCentres();
+    const tmp<pointField> tctrs(this->transform(mesh_.faceCentres()));
+    const pointField& ctrs = tctrs();
 
     forAll(ctrs, elemi)
     {
@@ -138,7 +139,7 @@ Foam::boxToFace::boxToFace
     const dictionary& dict
 )
 :
-    topoSetFaceSource(mesh),
+    topoSetFaceSource(mesh, dict),
     bbs_()
 {
     // Accept 'boxes', 'box' or 'min/max'

src/meshTools/topoSet/faceSources/boxToFace/boxToFace.H

@@ -50,7 +50,7 @@ Usage
         source      boxToFace;
 
         // Conditional mandatory entries
-        // Select either of the below
+        // Select one of the below
 
         // Option-1
         boxes

src/meshTools/topoSet/faceSources/cellToFace/cellToFace.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -68,23 +68,14 @@ Foam::cellToFace::cellActionNames_
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::cellToFace::combine
+template<class Selector>
+void Foam::cellToFace::combineImpl
 (
     topoSet& set,
     const bool add,
-    const word& setName
+    const Selector& cellLabels
 ) const
 {
-    // Load the set
-    if (!exists(mesh_.time().path()/topoSet::localPath(mesh_, setName)))
-    {
-        SeriousError<< "Cannot load set "
-            << setName << endl;
-    }
-
-    cellSet loadedSet(mesh_, setName);
-    const labelHashSet& cellLabels = loadedSet;
-
     if (option_ == ALL)
     {
         // Add all faces from cell
@@ -212,6 +203,36 @@ void Foam::cellToFace::combine
 }
 
 
+void Foam::cellToFace::combine
+(
+    topoSet& set,
+    const bool add,
+    const word& setName
+) const
+{
+    if (isZone_)
+    {
+        const labelList& cellLabels = mesh_.cellZones()[setName];
+
+        combineImpl(set, add, cellLabels);
+    }
+    else
+    {
+        // Load the set
+        if (!exists(mesh_.time().path()/topoSet::localPath(mesh_, setName)))
+        {
+            SeriousError<< "Cannot load set "
+                << setName << endl;
+        }
+
+        cellSet loadedSet(mesh_, setName);
+        const labelHashSet& cellLabels = loadedSet;
+
+        combineImpl(set, add, cellLabels);
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::cellToFace::cellToFace
@@ -222,7 +243,8 @@ Foam::cellToFace::cellToFace
 )
 :
     topoSetFaceSource(mesh),
-    names_(one{}, setName),
+    names_(Foam::one{}, setName),
+    isZone_(false),
     option_(option)
 {}
 
@@ -233,17 +255,11 @@ Foam::cellToFace::cellToFace
     const dictionary& dict
 )
 :
-    topoSetFaceSource(mesh),
+    topoSetFaceSource(mesh, dict),
     names_(),
+    isZone_(topoSetSource::readNames(dict, names_)),
     option_(cellActionNames_.get("option", dict))
-{
-    // Look for 'sets' or 'set'
-    if (!dict.readIfPresent("sets", names_))
-    {
-        names_.resize(1);
-        dict.readEntry("set", names_.front());
-    }
-}
+{}
 
 
 Foam::cellToFace::cellToFace
@@ -253,7 +269,8 @@ Foam::cellToFace::cellToFace
 )
 :
     topoSetFaceSource(mesh),
-    names_(one{}, word(checkIs(is))),
+    names_(Foam::one{}, word(checkIs(is))),
+    isZone_(false),
     option_(cellActionNames_.read(checkIs(is)))
 {}
 
@@ -270,7 +287,8 @@ void Foam::cellToFace::applyToSet
     {
         if (verbose_)
         {
-            Info<< "    Adding faces according to cell sets: "
+            Info<< "    Adding faces according to cell "
+                << (isZone_ ? "zones: " : "sets: ")
                 << flatOutput(names_) << nl;
         }
 
@@ -283,7 +301,8 @@ void Foam::cellToFace::applyToSet
     {
         if (verbose_)
         {
-            Info<< "    Removing faces according to cell sets: "
+            Info<< "    Removing faces according to cell "
+                << (isZone_ ? "zones: " : "sets: ")
                 << flatOutput(names_) << nl;
         }
 

src/meshTools/topoSet/faceSources/cellToFace/cellToFace.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -51,7 +51,7 @@ Usage
         option      <option>;
 
         // Conditional mandatory entries
-        // Select either of the below
+        // Select one of the below
 
         // Option-1
         sets
@@ -62,7 +62,18 @@ Usage
         );
 
         // Option-2
-        set <cellSetName>;
+        zones
+        (
+            <cellZoneName0>
+            <cellZoneName1>
+            ...
+        );
+
+        // Option-3
+        set     <cellSetName>;
+
+        // Option-4
+        zone    <cellZoneName>;
     }
     \endverbatim
 
@@ -90,16 +101,16 @@ Usage
       outside | Faces with only one neighbour in the cellSet
     \endverbatim
 
-    Options for the conditional mandatory entries:
+    Options for the conditional mandatory entries (in order of precedence):
     \verbatim
       Entry    | Description                    | Type     | Reqd   | Deflt
       sets     | Names of input cellSets        | wordList | choice | -
+      zones    | Names of input cellZones       | wordList | cond'l | -
       set      | Name of input cellSet          | word     | choice | -
+      zone     | Name of input cellZone         | word     | cond'l | -
     \endverbatim
 
 Note
-  - The order of precedence among the conditional mandatory entries from the
-    highest to the lowest is \c sets, and \c set.
   - The \c outside option applies to the cellSets individually.
 
 See also
@@ -149,9 +160,12 @@ private:
 
         static const Enum<cellAction> cellActionNames_;
 
-        //- Names of cellSets to use
+        //- Names of sets or zones to use
         wordList names_;
 
+        //- Is name a set or a zone
+        const bool isZone_;
+
         //- Selection type
         cellAction option_;
 
@@ -159,6 +173,15 @@ private:
     // Private Member Functions
 
         //- Depending on face to cell option add to or delete from cellSet.
+        template<class Selector>
+        void combineImpl
+        (
+            topoSet& set,
+            const bool add,
+            const Selector& cellLabels
+        ) const;
+
+        //- Depending on face to cell option add to or delete from set.
         void combine(topoSet& set, const bool add, const word& setName) const;
 
 

src/meshTools/topoSet/faceSources/clipPlaneToFace/clipPlaneToFace.C

@@ -69,7 +69,8 @@ void Foam::clipPlaneToFace::combine(topoSet& set, const bool add) const
 {
     // Face centres above the plane
 
-    const pointField& ctrs = mesh_.faceCentres();
+    const tmp<pointField> tctrs(this->transform(mesh_.faceCentres()));
+    const pointField& ctrs = tctrs();
 
     forAll(ctrs, elemi)
     {
@@ -102,12 +103,9 @@ Foam::clipPlaneToFace::clipPlaneToFace
     const dictionary& dict
 )
 :
-    clipPlaneToFace
-    (
-        mesh,
-        dict.get<vector>("point"),
-        dict.get<vector>("normal")
-    )
+    topoSetFaceSource(mesh, dict),
+    point_(dict.get<vector>("point")),
+    normal_(dict.get<vector>("normal"))
 {}
 
 

src/meshTools/topoSet/faceSources/cylinderToFace/cylinderToFace.C

@@ -68,7 +68,8 @@ Foam::topoSetSource::addToUsageTable Foam::cylinderToFace::usage_
 
 void Foam::cylinderToFace::combine(topoSet& set, const bool add) const
 {
-    const pointField& ctrs = mesh_.faceCentres();
+    const tmp<pointField> tctrs(this->transform(mesh_.faceCentres()));
+    const pointField& ctrs = tctrs();
 
     const vector axis = (point2_ - point1_);
     const scalar magAxis2 = magSqr(axis);
@@ -127,12 +128,12 @@ Foam::cylinderToFace::cylinderToFace
     const dictionary& dict
 )
 :
-    cylinderToFace
+    topoSetFaceSource(mesh, dict),
+    point1_(dict.getCompat<point>("point1", {{"p1", -2112}})),
+    point2_(dict.getCompat<point>("point2", {{"p2", -2112}})),
+    radius_(dict.getCompat<scalar>("radius", {{"outerRadius", -2112}})),
+    innerRadius_
     (
-        mesh,
-        dict.getCompat<point>("point1", {{"p1", -2112}}),
-        dict.getCompat<point>("point2", {{"p2", -2112}}),
-        dict.getCompat<scalar>("radius", {{"outerRadius", -2112}}),
         dict.getCheckOrDefault<scalar>("innerRadius", 0, scalarMinMax::ge(0))
     )
 {}

src/meshTools/topoSet/faceSources/faceToFace/faceToFace.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -60,7 +60,8 @@ Foam::faceToFace::faceToFace
 )
 :
     topoSetFaceSource(mesh),
-    names_(one{}, setName)
+    names_(Foam::one{}, setName),
+    isZone_(false)
 {}
 
 
@@ -70,16 +71,10 @@ Foam::faceToFace::faceToFace
     const dictionary& dict
 )
 :
-    topoSetFaceSource(mesh),
-    names_()
-{
-    // Look for 'sets' or 'set'
-    if (!dict.readIfPresent("sets", names_))
-    {
-        names_.resize(1);
-        dict.readEntry("set", names_.front());
-    }
-}
+    topoSetFaceSource(mesh, dict),
+    names_(),
+    isZone_(topoSetSource::readNames(dict, names_))
+{}
 
 
 Foam::faceToFace::faceToFace
@@ -89,7 +84,8 @@ Foam::faceToFace::faceToFace
 )
 :
     topoSetFaceSource(mesh),
-    names_(one{}, word(checkIs(is)))
+    names_(Foam::one{}, word(checkIs(is))),
+    isZone_(false)
 {}
 
 
@@ -105,30 +101,46 @@ void Foam::faceToFace::applyToSet
     {
         if (verbose_)
         {
-            Info<< "    Adding all elements of face sets: "
+            Info<< "    Adding all elements of face "
+                << (isZone_ ? "zones: " : "sets: ")
                 << flatOutput(names_) << nl;
         }
 
         for (const word& setName : names_)
         {
-            faceSet loadedSet(mesh_, setName);
+            if (isZone_)
+            {
+                set.addSet(mesh_.faceZones()[setName]);
+            }
+            else
+            {
+                faceSet loadedSet(mesh_, setName);
 
-            set.addSet(loadedSet);
+                set.addSet(loadedSet);
+            }
         }
     }
     else if (action == topoSetSource::SUBTRACT)
     {
         if (verbose_)
         {
-            Info<< "    Removing all elements of face sets: "
+            Info<< "    Removing all elements of face "
+                << (isZone_ ? "zones: " : "sets: ")
                 << flatOutput(names_) << nl;
         }
 
         for (const word& setName : names_)
         {
-            faceSet loadedSet(mesh_, setName);
+            if (isZone_)
+            {
+                set.subtractSet(mesh_.faceZones()[setName]);
+            }
+            else
+            {
+                faceSet loadedSet(mesh_, setName);
 
-            set.subtractSet(loadedSet);
+                set.subtractSet(loadedSet);
+            }
         }
     }
 }

src/meshTools/topoSet/faceSources/faceToFace/faceToFace.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -50,7 +50,7 @@ Usage
         source      faceToFace;
 
         // Conditional mandatory entries
-        // Select either of the below
+        // Select one of the below
 
         // Option-1
         sets
@@ -61,7 +61,18 @@ Usage
         );
 
         // Option-2
+        zones
+        (
+            <faceZoneName0>
+            <faceZoneName1>
+            ...
+        );
+
+        // Option-3
         set     <faceSetName>;
+
+        // Option-4
+        zone    <faceZoneName>;
     }
     \endverbatim
 
@@ -81,11 +92,13 @@ Usage
       subtract | Remove selected faces from this faceSet
     \endverbatim
 
-    Options for the conditional mandatory entries:
+    Options for the conditional mandatory entries (in order of precedence):
     \verbatim
       Entry    | Description                | Type     | Req'd  | Dflt
       sets     | Names of input faceSets    | wordList | cond'l | -
+      zones    | Names of input faceZones   | wordList | cond'l | -
       set      | Name of input faceSet      | word     | cond'l | -
+      zone     | Name of input faceZone     | word     | cond'l | -
     \endverbatim
 
 Note
@@ -101,8 +114,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef faceToFace_H
-#define faceToFace_H
+#ifndef Foam_faceToFace_H
+#define Foam_faceToFace_H
 
 #include "topoSetFaceSource.H"
 
@@ -124,9 +137,12 @@ class faceToFace
         //- Add usage string
         static addToUsageTable usage_;
 
-        //- Names of faceSets to use
+        //- Names of sets or zones to use
         wordList names_;
 
+        //- Is name a set or a zone
+        const bool isZone_;
+
 
 public:
 

src/meshTools/topoSet/faceSources/holeToFace/holeToFace.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020-2022 OpenCFD Ltd.
+    Copyright (C) 2020-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -1095,7 +1095,7 @@ Foam::holeToFace::holeToFace
     const dictionary& dict
 )
 :
-    topoSetFaceSource(mesh),
+    topoSetFaceSource(mesh, dict),
     zonePoints_(dict.get<List<pointField>>("points")),
     blockedFaceNames_(),
     blockedCellNames_(),
@@ -1149,7 +1149,7 @@ void Foam::holeToFace::applyToSet
     for (const word& setName : blockedFaceNames_)
     {
         const faceSet loadedSet(mesh_, setName);
-        isBlockedFace.set(loadedSet.toc());
+        isBlockedFace.setMany(loadedSet.begin(), loadedSet.end());
     }
 
     // Optional initial blocked cells
@@ -1159,7 +1159,7 @@ void Foam::holeToFace::applyToSet
         for (const word& setName : blockedCellNames_)
         {
             const cellSet loadedSet(mesh_, setName);
-            isCandidateCell.set(loadedSet.toc());
+            isCandidateCell.setMany(loadedSet.begin(), loadedSet.end());
         }
     }
     else

src/meshTools/topoSet/faceSources/labelToFace/labelToFace.C

@@ -94,11 +94,8 @@ Foam::labelToFace::labelToFace
     const dictionary& dict
 )
 :
-    labelToFace
-    (
-        mesh,
-        dict.get<labelList>("value")
-    )
+    topoSetFaceSource(mesh, dict),
+    labels_(dict.get<labelList>("value"))
 {}
 
 

src/meshTools/topoSet/faceSources/normalToFace/normalToFace.C

@@ -100,12 +100,9 @@ Foam::normalToFace::normalToFace
 
 Foam::normalToFace::normalToFace(const polyMesh& mesh, const dictionary& dict)
 :
-    normalToFace
-    (
-        mesh,
-        dict.get<vector>("normal"),
-        dict.get<scalar>("cos")
-    )
+    topoSetFaceSource(mesh, dict),
+    normal_(dict.get<vector>("normal")),
+    tol_(dict.get<scalar>("cos"))
 {
     setNormal();
 }

src/meshTools/topoSet/faceSources/patchToFace/patchToFace.C

@@ -121,7 +121,7 @@ Foam::patchToFace::patchToFace
     const dictionary& dict
 )
 :
-    topoSetFaceSource(mesh),
+    topoSetFaceSource(mesh, dict),
     selectedPatches_()
 {
     // Look for 'patches' and 'patch', but accept 'name' as well

src/meshTools/topoSet/faceSources/patchToFace/patchToFace.H

@@ -49,7 +49,7 @@ Usage
         source      patchToFace;
 
         // Conditional mandatory entries
-        // Select either of the below
+        // Select one of the below
 
         // Option-1
         patches

src/meshTools/topoSet/faceSources/pointToFace/pointToFace.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -82,18 +82,14 @@ Foam::pointToFace::pointActionNames_
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::pointToFace::combine
+template<class Selector>
+void Foam::pointToFace::combineImpl
 (
     topoSet& set,
     const bool add,
-    const word& setName
+    const Selector& pointLabels
 ) const
 {
-    // Load the set
-    pointSet loadedSet(mesh_, setName);
-
-    const labelHashSet& pointLabels = loadedSet;
-
     if (option_ == ANY)
     {
         // Add faces with any point in loadedSet
@@ -160,6 +156,31 @@ void Foam::pointToFace::combine
 }
 
 
+void Foam::pointToFace::combine
+(
+    topoSet& set,
+    const bool add,
+    const word& setName
+) const
+{
+    if (isZone_)
+    {
+        const labelList& pointLabels = mesh_.pointZones()[setName];
+
+        combineImpl(set, add, pointLabels);
+    }
+    else
+    {
+        // Load the set
+        pointSet loadedSet(mesh_, setName);
+
+        const labelHashSet& pointLabels = loadedSet;
+
+        combineImpl(set, add, pointLabels);
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::pointToFace::pointToFace
@@ -170,7 +191,8 @@ Foam::pointToFace::pointToFace
 )
 :
     topoSetFaceSource(mesh),
-    names_(one{}, setName),
+    names_(Foam::one{}, setName),
+    isZone_(false),
     option_(option)
 {}
 
@@ -181,17 +203,11 @@ Foam::pointToFace::pointToFace
     const dictionary& dict
 )
 :
-    topoSetFaceSource(mesh),
+    topoSetFaceSource(mesh, dict),
     names_(),
+    isZone_(topoSetSource::readNames(dict, names_)),
     option_(pointActionNames_.get("option", dict))
-{
-    // Look for 'sets' or 'set'
-    if (!dict.readIfPresent("sets", names_))
-    {
-        names_.resize(1);
-        dict.readEntry("set", names_.front());
-    }
-}
+{}
 
 
 Foam::pointToFace::pointToFace
@@ -201,7 +217,8 @@ Foam::pointToFace::pointToFace
 )
 :
     topoSetFaceSource(mesh),
-    names_(one{}, word(checkIs(is))),
+    names_(Foam::one{}, word(checkIs(is))),
+    isZone_(false),
     option_(pointActionNames_.read(checkIs(is)))
 {}
 
@@ -218,7 +235,8 @@ void Foam::pointToFace::applyToSet
     {
         if (verbose_)
         {
-            Info<< "    Adding faces according to point sets: "
+            Info<< "    Adding faces according to point "
+                << (isZone_ ? "zones: " : "sets: ")
                 << flatOutput(names_) << nl;
         }
 
@@ -231,7 +249,8 @@ void Foam::pointToFace::applyToSet
     {
         if (verbose_)
         {
-            Info<< "    Removing faces according to point sets: "
+            Info<< "    Removing faces according to point "
+                << (isZone_ ? "zones: " : "sets: ")
                 << flatOutput(names_) << nl;
         }