COMP: finiteArea: dangling reference warning

- can be useful for various orientation-related geometry or mesh
  manipulations during pre-/post-processing:
  * combine with linearDirection to achieve better extrusion results.
  * orientation of transformations, blockMesh, result projections, ...

STYLE: minor code modernizations

Co-authored-by: Mark Olesen <>

.gitlab/issue_templates/bug.md

@@ -49,10 +49,10 @@
 
 <!--
   Providing details of your set-up can help us identify any issues, e.g.
-  OpenFOAM version : v2312|v2306|v2212|v2206|v2112 etc
-  Operating system : ubuntu|openSUSE|centos etc
+  OpenFOAM version : v2406|v2312|v2306|v2212|v2206 etc
+  Operating system : ubuntu|openSUSE|RedHat etc
   Hardware info    : any info that may help?
-  Compiler         : gcc|intel|clang etc
+  Compiler         : gcc|clang etc
  -->
 
 - OpenFOAM version :

.gitmodules

@@ -1,8 +1,8 @@
 [submodule "cfmesh"]
-	path = modules/cfmesh
+	path = plugins/cfmesh
 	url = https://develop.openfoam.com/Community/integration-cfmesh.git
 [submodule "avalanche"]
-	path = modules/avalanche
+	path = plugins/avalanche
 	url = https://develop.openfoam.com/Community/avalanche.git
 [submodule "adios"]
 	path = modules/adios
@@ -18,5 +18,8 @@
 	path = modules/external-solver
 	url = https://develop.openfoam.com/Modules/external-solver.git
 [submodule "turbulence-community"]
-	path = modules/turbulence-community
+	path = plugins/turbulence-community
 	url = https://gitlab.com/openfoam/community/tc-turbulence/turbulence-community.git
+[submodule "plugins/data-community"]
+	path = plugins/data-community
+	url = https://gitlab.com/openfoam/community/sig-data-modelling/data-community.git

META-INFO/api-info

@@ -1,2 +1,2 @@
-api=2402
-patch=240522
+api=2406
+patch=0

README.md

@@ -40,9 +40,9 @@ Violations of the Trademark are monitored, and will be duly prosecuted.
 
 If OpenFOAM has already been compiled on your system, simply source
 the appropriate `etc/bashrc` or `etc/cshrc` file and get started.
-For example, for the OpenFOAM-v2312 version:
+For example, for the OpenFOAM-v2406 version:
 ```
-source /installation/path/OpenFOAM-v2312/etc/bashrc
+source /installation/path/OpenFOAM-v2406/etc/bashrc
 ```
 
 ## Compiling OpenFOAM
@@ -127,8 +127,8 @@ These 3rd-party sources are normally located in a directory parallel
 to the OpenFOAM directory. For example,
 ```
 /path/parent
-|-- OpenFOAM-v2312
-\-- ThirdParty-v2312
+|-- OpenFOAM-v2406
+\-- ThirdParty-v2406
 ```
 There are, however, many cases where this simple convention is inadequate:
 
@@ -136,7 +136,7 @@ There are, however, many cases where this simple convention is inadequate:
   operating system or cluster installation provides it)
 
 * When we have changed the OpenFOAM directory name to some arbitrary
-  directory name, e.g. openfoam-sandbox2312, etc..
+  directory name, e.g. openfoam-sandbox2406, etc..
 
 * When we would like any additional 3rd party software to be located
   inside of the OpenFOAM directory to ensure that the installation is
@@ -156,9 +156,9 @@ when locating the ThirdParty directory with the following precedence:
 2. PREFIX/ThirdParty-VERSION
    * this corresponds to the traditional approach
 3. PREFIX/ThirdParty-vAPI
-   * allows for an updated value of VERSION, *eg*, `v2312-myCustom`,
+   * allows for an updated value of VERSION, *eg*, `v2406-myCustom`,
      without requiring a renamed ThirdParty. The API value would still
-     be `2312` and the original `ThirdParty-v2312/` would be found.
+     be `2406` and the original `ThirdParty-v2406/` would be found.
 4. PREFIX/ThirdParty-API
    * same as the previous example, but using an unadorned API value.
 5. PREFIX/ThirdParty-common
@@ -213,4 +213,4 @@ ThirdParty directory will contain either an `Allwmake` file or a
 - [Governance](http://www.openfoam.com/governance/), [Governance Projects](https://www.openfoam.com/governance/projects)
 - [Contacting OpenCFD](http://www.openfoam.com/contact/)
 
-Copyright 2016-2023 OpenCFD Ltd
+Copyright 2016-2024 OpenCFD Ltd

applications/solvers/combustion/reactingFoam/EEqn.H

@@ -21,6 +21,11 @@
       + fvOptions(rho, he)
     );
 
+    if (MRF.active())
+    {
+        EEqn += fvc::div(MRF.phi(), p);
+    }
+
     EEqn.relax();
 
     fvOptions.constrain(EEqn);

applications/solvers/compressible/rhoPimpleFoam/EEqn.H

@@ -20,6 +20,11 @@
         fvOptions(rho, he)
     );
 
+    if (MRF.active())
+    {
+        EEqn += fvc::div(MRF.phi(), p);
+    }
+
     EEqn.relax();
 
     fvOptions.constrain(EEqn);

applications/solvers/compressible/rhoSimpleFoam/EEqn.H

@@ -14,6 +14,11 @@
         fvOptions(rho, he)
     );
 
+    if (MRF.active())
+    {
+        EEqn += fvc::div(MRF.phi(), p);
+    }
+
     EEqn.relax();
 
     fvOptions.constrain(EEqn);

applications/solvers/compressible/sonicFoam/EEqn.H

@@ -9,6 +9,11 @@
         fvOptions(rho, e)
     );
 
+    if (MRF.active())
+    {
+        EEqn += fvc::div(MRF.phi(), p);
+    }
+
     EEqn.relax();
 
     fvOptions.constrain(EEqn);

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/icoReactingMultiphaseInterFoam.C

@@ -130,10 +130,17 @@ int main(int argc, char *argv[])
                     }
                 }
             }
-            #include "UEqn.H"
+
             #include "YEqns.H"
             #include "TEqn.H"
 
+            if (pimple.frozenFlow())
+            {
+                continue;
+            }
+
+            #include "UEqn.H"
+
             // --- Pressure corrector loop
             while (pimple.correct())
             {

applications/test/exprValue1/Test-exprValue1.cxx

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021-2023 OpenCFD Ltd.
+    Copyright (C) 2021-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -21,6 +21,8 @@ Description
 #include "argList.H"
 #include "IOstreams.H"
 #include "ITstream.H"
+#include "OTstream.H"
+#include "SpanStream.H"
 #include "exprValue.H"
 #include "Pstream.H"
 
@@ -34,6 +36,41 @@ void printInfo(const expressions::exprValue& val)
 }
 
 
+void write_read(const expressions::exprValue& val)
+{
+    OCharStream os;
+    os << val;
+
+    ISpanStream is(os.view());
+    expressions::exprValue val2;
+    is >> val2;
+
+    Pout<< "wrote " << os.count() << " chars: " << os.str() << nl;
+
+    printInfo(val);
+    printInfo(val2);
+    Pout<< "====" << nl;
+}
+
+
+tokenList tokens_of(const expressions::exprValue& val)
+{
+    OTstream toks;
+    toks << val;
+
+    Pout<< "val with tokens: ";
+    toks.writeList(Pout, 0) << nl;
+
+    for (const auto& t : toks)
+    {
+        Pout<< "  " << t.info() << nl;
+    }
+    Pout<< nl;
+
+    return toks;
+}
+
+
 expressions::exprValue tryParse(const std::string& str)
 {
     expressions::exprValue val, val2;
@@ -90,6 +127,7 @@ int main(int argc, char *argv[])
 
     {
         expressions::exprValue value;
+        tokenList toks;
 
         Info<< "exprValue"
             << " sizeof:" << value.size_bytes()
@@ -99,21 +137,31 @@ int main(int argc, char *argv[])
 
         // Nothing
         printInfo(value);
+        toks = tokens_of(value);
+        write_read(value);
 
         value.set(scalar(100));
-        printInfo(value);
+        printInfo(value); write_read(value); toks = tokens_of(value);
+
+        value.set(scalar(100.01));
+        printInfo(value); write_read(value); toks = tokens_of(value);
 
         value.set(vector(1,2,3));
-        printInfo(value);
+        printInfo(value); write_read(value); toks = tokens_of(value);
 
         value = vector(4,5,6);
-        printInfo(value);
+        printInfo(value); write_read(value); toks = tokens_of(value);
 
         value = Zero;
-        printInfo(value);
+        printInfo(value); write_read(value); toks = tokens_of(value);
 
         value.clear();
+        printInfo(value); write_read(value); toks = tokens_of(value);
+
+        value.set<bool>(true);
+
         printInfo(value);
+        printInfo(value); write_read(value); toks = tokens_of(value);
 
         if (UPstream::parRun())
         {

applications/test/exprValue2/Make/files

@@ -1,5 +1,3 @@
-exprValueFieldTag.cxx
-
 Test-exprValue2.cxx
 
 EXE = $(FOAM_USER_APPBIN)/Test-exprValue2

applications/test/exprValue2/Test-exprValue2.cxx

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2023 OpenCFD Ltd.
+    Copyright (C) 2023-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -36,7 +36,6 @@ Description
 #include "vectorField.H"
 #include "DynamicList.H"
 #include "Random.H"
-#include "exprValue.H"
 #include "exprValueFieldTag.H"
 
 using namespace Foam;
@@ -61,27 +60,81 @@ int main(int argc, char *argv[])
 
     #include "setRootCase.H"
 
+    DynamicList<fieldTag> allTags;
+
     {
         scalarField fld1(20);
         scalarField fld2a(20, Zero);
         scalarField fld2b(10, 3.10);
         scalarField fld3;
 
-        forAll(fld1, i)
+        for (auto& val : fld1)
         {
-            fld1[i] = rnd.position<scalar>(0, 20);
+            val = rnd.position<scalar>(0, 20);
+        }
+
+        if (!UPstream::master())
+        {
+            fld2b.resize(5);
+            fld2b *= 2;
         }
 
         fieldTag tag1(fld1.begin(), fld1.end());
         fieldTag tag2a(fld2a.begin(), fld2a.end());
         fieldTag tag2b(fld2b.begin(), fld2b.end());
         fieldTag tag3(fld3.begin(), fld3.end());
+        fieldTag tag4(fld3.begin(), fld3.end());
 
         printInfo(tag1) << nl;
         printInfo(tag2a) << nl;
         printInfo(tag2b) << nl;
         printInfo(tag3) << nl;
 
+        {
+            Pout<< "Test reduce" << nl;
+
+            fieldTag work(fld2b.begin(), fld2b.end());
+
+            Pout<< "Before" << nl;
+            printInfo(work) << nl;
+
+            work.reduce();
+
+            Pout<< "After" << nl;
+            printInfo(work) << nl;
+            Pout<< "====" << nl;
+        }
+
+        allTags.clear();
+        allTags.push_back(tag1);
+        allTags.push_back(tag2a);
+        allTags.push_back(tag2b);
+        allTags.push_back(tag3);
+        allTags.push_back(tag4);
+        allTags.push_back(fieldTag::make_empty<tensor>());
+
+
+        // Add some other types
+        {
+            vectorField vfld2a(20, vector::uniform(1.23));
+
+            allTags.emplace_back
+            (
+                vfld2a.begin(),
+                vfld2a.end()
+            );
+            allTags.emplace_back(vector(1.01, 2.02, 3.03));
+            allTags.emplace_back(12.4);
+
+            allTags.emplace_back().set_value(vector::uniform(2.0));
+            allTags.back().set_empty();
+        }
+        Info<< "all tags: " << allTags << nl;
+
+        Foam::sort(allTags);
+
+        Info<< "sorted: " << allTags << nl;
+
         fieldTag result;
 
         result = fieldTag::combineOp{}(tag1, tag2a);

applications/test/exprValue2/exprValueFieldTag.cxx

@@ -1,160 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2023 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "exprValueFieldTag.H"
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-bool Foam::expressions::exprValueFieldTag::empty() const noexcept
-{
-    return
-    (
-        uniformity_ == Foam::Detail::ListPolicy::uniformity::EMPTY
-    );
-}
-
-
-bool Foam::expressions::exprValueFieldTag::is_uniform() const noexcept
-{
-    return
-    (
-        uniformity_ == Foam::Detail::ListPolicy::uniformity::UNIFORM
-    );
-}
-
-
-bool Foam::expressions::exprValueFieldTag::is_nonuniform() const noexcept
-{
-    return
-    (
-        uniformity_ == Foam::Detail::ListPolicy::uniformity::NONUNIFORM
-    );
-}
-
-
-bool Foam::expressions::exprValueFieldTag::equal
-(
-    const exprValueFieldTag& rhs
-) const
-{
-    return (value_ == rhs.value_);
-}
-
-
-void Foam::expressions::exprValueFieldTag::set_nouniform() noexcept
-{
-    uniformity_ = Foam::Detail::ListPolicy::uniformity::NONUNIFORM;
-    value_ = Foam::zero{};
-}
-
-
-void Foam::expressions::exprValueFieldTag::combine
-(
-    const exprValueFieldTag& b
-)
-{
-    if (b.empty())
-    {
-        // no-op
-        return;
-    }
-
-    exprValueFieldTag& a = *this;
-
-    if (a.empty())
-    {
-        a = b;
-    }
-    else if (a.is_nonuniform())
-    {
-        // Already non-uniform/mixed
-        // a.uniformity_ |= b.uniformity_;
-
-        a.value_ = Foam::zero{};
-    }
-    else if (a.is_uniform() && b.is_uniform())
-    {
-        // Both are uniform, but are they the same value?
-        if (!a.equal(b))
-        {
-            a.set_nouniform();
-        }
-    }
-}
-
-
-// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
-
-void Foam::expressions::exprValueFieldTag::read(Istream& is)
-{
-    label uniformTag;
-
-    is >> uniformTag;
-    uniformity_ = int(uniformTag);
-    value_.read(is);
-}
-
-
-void Foam::expressions::exprValueFieldTag::write(Ostream& os) const
-{
-    os << label(uniformity_);
-    value_.write(os, false);  // No pruning
-}
-
-
-void Foam::expressions::exprValueFieldTag::print(Ostream& os) const
-{
-    os  << "{ uniform:"
-        << label(uniformity_)
-        << " type:" << label(value_.typeCode())
-        << " value: " << value_ << " }";
-}
-
-
-Foam::Istream& Foam::operator>>
-(
-    Istream& is,
-    expressions::exprValueFieldTag& tag
-)
-{
-    tag.read(is);
-    return is;
-}
-
-
-Foam::Ostream& Foam::operator<<
-(
-    Ostream& os,
-    const expressions::exprValueFieldTag& tag
-)
-{
-    tag.write(os);
-    return os;
-}
-
-
-// ************************************************************************* //

applications/utilities/mesh/manipulation/transformPoints/transformPoints.C

@@ -7,6 +7,7 @@
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
     Copyright (C) 2017-2024 OpenCFD Ltd.
+    Copyright (C) 2024 Haakan Nilsson
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -71,6 +72,9 @@ Usage
     The any or all of the three options may be specified and are processed
     in the above order.
 
+    -cylToCart (originVector axisVector directionVector)
+        Tranform cylindrical coordinates to cartesian coordinates
+
     With -rotateFields (in combination with -rotate/yawPitchRoll/rollPitchYaw)
     it will also read & transform vector & tensor fields.
 
@@ -93,6 +97,7 @@ Note
 #include "transformGeometricField.H"
 #include "axisAngleRotation.H"
 #include "EulerCoordinateRotation.H"
+#include "cylindricalCS.H"
 
 using namespace Foam;
 using namespace Foam::coordinateRotations;
@@ -326,6 +331,13 @@ int main(int argc, char *argv[])
         "Scale by the specified amount - Eg, for uniform [mm] to [m] scaling "
         "use either '(0.001 0.001 0.001)' or simply '0.001'"
     );
+    argList::addOption
+    (
+        "cylToCart",
+        "(originVec axisVec directionVec)",
+        "Tranform cylindrical coordinates to cartesian coordinates"
+    );
+
 
     // Compatibility with surfaceTransformPoints
     argList::addOptionCompat("scale", {"write-scale", 0});
@@ -348,7 +360,8 @@ int main(int argc, char *argv[])
             "rotate-z",
             "rollPitchYaw",
             "yawPitchRoll",
-            "scale"
+            "scale",
+            "cylToCart"
         });
 
         if (!args.count(operationNames))
@@ -552,6 +565,23 @@ int main(int argc, char *argv[])
         // Output scaling
         applyScaling(points, getScalingOpt("scale", args));
 
+        if (args.found("cylToCart"))
+        {
+            vectorField n1n2(args.lookup("cylToCart")());
+            n1n2[1].normalise();
+            n1n2[2].normalise();
+
+            cylindricalCS ccs
+                (
+                    "ccs",
+                    n1n2[0],
+                    n1n2[1],
+                    n1n2[2]
+                );
+
+            points = ccs.globalPosition(points);
+        }
+
         // More precision (for points data)
         IOstream::minPrecision(10);
 

applications/utilities/surface/surfaceInertia/surfaceInertia.C

@@ -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.
@@ -161,25 +161,22 @@ int main(int argc, char *argv[])
         // rather than tensors to allow indexed permutation.
 
         // Cartesian basis vectors - right handed orthogonal triplet
-        List<vector> cartesian(3);
+        FixedList<vector, 3> cartesian;
 
         cartesian[0] = vector(1, 0, 0);
         cartesian[1] = vector(0, 1, 0);
         cartesian[2] = vector(0, 0, 1);
 
-        // Principal axis basis vectors - right handed orthogonal
-        // triplet
-        List<vector> principal(3);
+        // Principal axis basis vectors - right handed orthogonal triplet
+        FixedList<vector, 3> principal;
 
         principal[0] = eVec.x();
         principal[1] = eVec.y();
         principal[2] = eVec.z();
 
-        scalar maxMagDotProduct = -GREAT;
-
         // Matching axis indices, first: cartesian, second:principal
-
         Pair<label> match(-1, -1);
+        scalar maxMagDotProduct = -GREAT;
 
         forAll(cartesian, cI)
         {
@@ -187,12 +184,11 @@ int main(int argc, char *argv[])
             {
                 scalar magDotProduct = mag(cartesian[cI] & principal[pI]);
 
-                if (magDotProduct > maxMagDotProduct)
+                if (maxMagDotProduct < magDotProduct)
                 {
                     maxMagDotProduct = magDotProduct;
 
                     match.first() = cI;
-
                     match.second() = pI;
                 }
             }
@@ -208,7 +204,7 @@ int main(int argc, char *argv[])
             // Invert the best match direction and swap the order of
             // the other two vectors
 
-            List<vector> tPrincipal = principal;
+            FixedList<vector, 3> tPrincipal = principal;
 
             tPrincipal[match.second()] *= -1;
 
@@ -237,7 +233,7 @@ int main(int argc, char *argv[])
 
             permutationDelta += 3;
 
-            List<vector> tPrincipal = principal;
+            FixedList<vector, 3> tPrincipal = principal;
 
             vector tEVal = eVal;
 
@@ -253,10 +249,9 @@ int main(int argc, char *argv[])
             eVal = tEVal;
         }
 
-        label matchedAlready = match.first();
-
-        match =Pair<label>(-1, -1);
+        const label matchedAlready = match.first();
 
+        match = Pair<label>(-1, -1);
         maxMagDotProduct = -GREAT;
 
         forAll(cartesian, cI)
@@ -275,12 +270,11 @@ int main(int argc, char *argv[])
 
                 scalar magDotProduct = mag(cartesian[cI] & principal[pI]);
 
-                if (magDotProduct > maxMagDotProduct)
+                if (maxMagDotProduct < magDotProduct)
                 {
                     maxMagDotProduct = magDotProduct;
 
                     match.first() = cI;
-
                     match.second() = pI;
                 }
             }
@@ -295,7 +289,7 @@ int main(int argc, char *argv[])
         {
             principal[match.second()] *= -1;
 
-            List<vector> tPrincipal = principal;
+            FixedList<vector, 3> tPrincipal = principal;
 
             tPrincipal[(matchedAlready + 1) % 3] =
                 principal[(matchedAlready + 2) % 3]*-sense;
@@ -335,15 +329,16 @@ int main(int argc, char *argv[])
         showTransform = false;
     }
 
-    // calculate the total surface area
 
+    // Calculate total surface area and average normal vector
     scalar surfaceArea = 0;
+    vector averageNormal(Zero);
 
     forAll(surf, facei)
     {
         const labelledTri& f = surf[facei];
 
-        if (f[0] == f[1] || f[0] == f[2] || f[1] == f[2])
+        if (!f.good())
         {
             WarningInFunction
                << "Illegal triangle " << facei << " vertices " << f
@@ -351,20 +346,23 @@ int main(int argc, char *argv[])
         }
         else
         {
-            surfaceArea += triPointRef
-            (
-                surf.points()[f[0]],
-                surf.points()[f[1]],
-                surf.points()[f[2]]
-            ).mag();
+            triPointRef tri = f.tri(surf.points());
+
+            surfaceArea += tri.mag();
+            averageNormal += tri.areaNormal();
         }
     }
 
+    // The unit normal (area-averaged)
+    averageNormal.normalise();
+
+
     Info<< nl << setprecision(12)
         << "Density: " << density << nl
         << "Mass: " << m << nl
         << "Centre of mass: " << cM << nl
         << "Surface area: " << surfaceArea << nl
+        << "Average normal: " << averageNormal << nl
         << "Inertia tensor around centre of mass: " << nl << J << nl
         << "eigenValues (principal moments): " << eVal << nl
         << "eigenVectors (principal axes): " << nl
@@ -398,21 +396,26 @@ int main(int argc, char *argv[])
             << endl;
     }
 
-    OFstream str("axes.obj");
 
-    Info<< nl << "Writing scaled principal axes at centre of mass of "
-        << surfFileName <<  " to " << str.name() << endl;
-
-    scalar scale = mag(cM - surf.points()[0])/eVal.component(findMin(eVal));
-
-    meshTools::writeOBJ(str, cM);
-    meshTools::writeOBJ(str, cM + scale*eVal.x()*eVec.x());
-    meshTools::writeOBJ(str, cM + scale*eVal.y()*eVec.y());
-    meshTools::writeOBJ(str, cM + scale*eVal.z()*eVec.z());
-
-    for (label i = 1; i < 4; i++)
+    // Write (scaled) principal axes at centre of mass
     {
-        str << "l " << 1 << ' ' << i + 1 << endl;
+        OFstream str("axes.obj");
+
+        Info<< nl << "Writing scaled principal axes at centre of mass of "
+            << surfFileName <<  " to " << str.name() << endl;
+
+        scalar scale =
+            mag(cM - surf.points()[0])/eVal.component(findMin(eVal));
+
+        meshTools::writeOBJ(str, cM);
+        meshTools::writeOBJ(str, cM + scale*eVal.x()*eVec.x());
+        meshTools::writeOBJ(str, cM + scale*eVal.y()*eVec.y());
+        meshTools::writeOBJ(str, cM + scale*eVal.z()*eVec.z());
+
+        for (label i = 1; i < 4; i++)
+        {
+            str << "l " << 1 << ' ' << i + 1 << endl;
+        }
     }
 
     Info<< "\nEnd\n" << endl;

applications/utilities/surface/surfaceTransformPoints/surfaceTransformPoints.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2017-2022 OpenCFD Ltd.
+    Copyright (C) 2017-2022,2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -58,6 +58,7 @@ Description
 #include "axisAngleRotation.H"
 #include "EulerCoordinateRotation.H"
 #include "MeshedSurfaces.H"
+#include "cylindricalCS.H"
 
 using namespace Foam;
 using namespace Foam::coordinateRotations;
@@ -257,6 +258,12 @@ int main(int argc, char *argv[])
         "type",
         "Output format (default: use file extension)"
     );
+    argList::addOption
+    (
+        "cylToCart",
+        "(originVec axisVec directionVec)",
+        "Tranform cylindrical coordinates to cartesian coordinates"
+    );
 
     // Backward compatibility and with transformPoints
     argList::addOptionCompat("write-scale", {"scale", -2006});
@@ -277,7 +284,8 @@ int main(int argc, char *argv[])
             "rollPitchYaw",
             "yawPitchRoll",
             "read-scale",
-            "write-scale"
+            "write-scale",
+            "cylToCart"
         });
 
         if (!args.count(operationNames))
@@ -461,6 +469,24 @@ int main(int argc, char *argv[])
     // Output scaling
     applyScaling(points, getScalingOpt("write-scale", args));
 
+    // Conversion to cylindrical coords
+    if (args.found("cylToCart"))
+    {
+        vectorField n1n2(args.lookup("cylToCart")());
+        n1n2[1].normalise();
+        n1n2[2].normalise();
+
+        cylindricalCS ccs
+            (
+                "ccs",
+                n1n2[0],
+                n1n2[1],
+                n1n2[2]
+            );
+
+        points = ccs.globalPosition(points);
+    }
+
     surf1.movePoints(points);
     surf1.write(exportName, writeFileType);
 

bin/tools/foamPackRelease

@@ -6,7 +6,7 @@
 #   \\  /    A nd           | www.openfoam.com
 #    \\/     M anipulation  |
 #------------------------------------------------------------------------------
-#     Copyright (C) 2019-2023 OpenCFD Ltd.
+#     Copyright (C) 2019-2024 OpenCFD Ltd.
 #------------------------------------------------------------------------------
 # License
 #     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -48,8 +48,16 @@ options:
   -name=NAME        Stem for tar-file (default: auto)
   -output=DIR       Output directory (default: ".")
   -prefix=NAME      Prefix directory within tar-file (default: auto)
-  -pkg-modules      Only package submodules - exclude OpenFOAM
-  -no-modules       Exclude submodules
+  -pkg-modules      Package 'modules' exclusively (no OpenFOAM)
+  -pkg-plugins      Package 'plugins' exclusively (no OpenFOAM)
+  -no-extras        Exclude 'modules, plugins,...' from source pack
+  -no-modules       Exclude 'modules' from source pack (default: off)
+  -no-plugins       Exclude 'plugins' from source pack (default: on)
+  -all-extras       Include 'modules, plugins,...' into source pack
+  -with-modules     Include 'modules' into source pack (default: on)
+  -with-plugins     Include 'plugins' into source pack (default: off)
+  -modules=name1,.. Include specifed 'modules' into source pack
+  -plugins=name1,.. Include specifed 'plugins' into source pack
   -no-patch         Ignore '_patch' number for output tar-file
   -no-prefix        Do not prefix subdirectory
   -no-compress      Disable compression
@@ -67,7 +75,7 @@ Script generator for packing OpenFOAM sources and submodules.
 Eg,
 
     $Script -output=some-dir origin/master > create-tar-file
-    sh ./create-tar-file
+    bash ./create-tar-file
 
     $Script -tgz origin/master | bash
 
@@ -90,11 +98,16 @@ die()
 
 
 #-------------------------------------------------------------------------------
-outputDir="."
+outputDir='.'
 versionSeparator='_'
 withPatchNum=true
-unset compress packageApi withSource withModules prefixDir tarName
-unset gitbase
+# Default selections
+select_source=true
+select_modules=true
+select_plugins=false
+
+unset compress gitbase packageApi prefixDir tarName
+
 
 # Cleanup tarName to remove trailing '.tar', detect compression etc
 cleanTarName() {
@@ -117,7 +130,7 @@ cleanTarName() {
 while [ "$#" -gt 0 ]
 do
     case "$1" in
-    -h | -help*)
+    -h | -help* | --help*)
         printHelp
         ;;
     -debian | -debian=[0-9]*)
@@ -146,13 +159,43 @@ do
         prefixDir="${1#*=}"
         prefixDir="${prefixDir%/}"
         ;;
-    -pkg-modules)
-        withModules=true
-        withSource=false
+    -all-extras)
+        select_modules=true
+        select_plugins=true
+        ;;
+    -with-modules)
+        select_modules=true
+        ;;
+    -modules=*)
+        select_modules="${1#*=}"
+        : "${select_modules:=true}"
+        ;;
+    -with-plugins)
+        select_plugins=true
+        ;;
+    -plugins=*)
+        select_plugins="${1#*=}"
+        : "${select_plugins:=true}"
+        ;;
+    -no-extras)
+        select_modules=false
+        select_plugins=false
         ;;
     -no-modules)
-        withModules=false
-        withSource=true
+        select_modules=false
+        ;;
+    -no-plugins)
+        select_plugins=false
+        ;;
+    -pkg-modules)  # Package modules exclusively
+        select_modules=true
+        select_plugins=false
+        select_source=false
+        ;;
+    -pkg-plugins)  # Package plugins exclusively
+        select_modules=false
+        select_plugins=true
+        select_source=false
         ;;
     -no-patch)
         withPatchNum=false
@@ -291,9 +334,15 @@ fi
 if [ -z "$prefixDir" ]
 then
     prefixDir="OpenFOAM-v${packageApi}"
-    if [ "$withSource" = false ]
+    if [ "$select_source" = false ]
     then
-        prefixDir="OpenFOAM-modules-v${packageApi}"
+        # Either -pkg-modules or -pkg-plugins, not both
+        if [ "$select_modules" != false ]
+        then
+            prefixDir="OpenFOAM-modules-v${packageApi}"
+        else
+            prefixDir="OpenFOAM-plugins-v${packageApi}"
+        fi
     fi
 elif [ "$prefixDir" = false ]
 then
@@ -320,9 +369,15 @@ case "$tarName" in
 
 ('')
     tarName="OpenFOAM-v${packageApi}"
-    if [ "$withSource" = false ]
+    if [ "$select_source" = false ]
     then
-        tarName="OpenFOAM-modules-v${packageApi}"
+        # Either -pkg-modules or -pkg-plugins, not both
+        if [ "$select_modules" != false ]
+        then
+            tarName="OpenFOAM-modules-v${packageApi}"
+        else
+            tarName="OpenFOAM-plugins-v${packageApi}"
+        fi
     fi
 
     if [ "$withPatchNum" = false ]
@@ -362,48 +417,56 @@ echo 'tar -cf "$outputDir/$tarName.tar" -T /dev/null'
 echo
 echo 'buildInfo="${prefixDir}${prefixDir:+@}META-INFO@build-info"'
 echo 'manifest0="${prefixDir}${prefixDir:+@}META-INFO@manifest.txt"'
-echo 'manifest1="${prefixDir}${prefixDir:+@}META-INFO@modules-manifest.txt"'
-echo '#--------'
+echo 'manifest1="${prefixDir}${prefixDir:+@}META-INFO@manifest-modules.txt"'
+echo 'manifest2="${prefixDir}${prefixDir:+@}META-INFO@manifest-plugins.txt"'
 
 #------------------------------------------------------------------------------
-# Sort out particulars related to modules
-if [ "$withModules" = false ]
+# Sort out particulars related to modules, source
+
+if [ "$select_source" = false ]
 then
-    echo '# No modules'
-    echo 'unset manifest1'
+    echo 'unset buildInfo manifest0  # No source'
 fi
 
-if [ "$withSource" = false ]
+if [ "$select_modules" = false ]
 then
-    echo '# No OpenFOAM source (package modules exclusively)'
-    echo 'unset buildInfo'
-    echo 'unset manifest0'
+    echo 'unset manifest1  # No modules'
 fi
 
+if [ "$select_plugins" = false ]
+then
+    echo 'unset manifest2  # No plugins'
+fi
+
+echo '#--------'
 echo 'set -x'
 echo
 
 #------------------------------------------------------------------------------
-# OpenFOAM sources
-if [ "$withSource" != false ]
+# OpenFOAM sources (unless explicitly excluded)
+
+if [ "$select_source" != false ]
 then
     echo 'git -c tar.umask=0022 archive --format=tar ${prefixDir:+--prefix="$prefixDir/"} -o "$outputDir/$tarName.tar" "$head"'
 
     # Tag build information with underscore to distinguish from "real" build
     # information when git is available.
-    echo 'echo build="${build:+_}$build" > "$outputDir/$buildInfo"'
-
+    echo echo 'build="${build:+_}$build" > "$outputDir/$buildInfo"'
+    echo
+    echo '# source'
+    echo 'manifestFile="$manifest0"'
     echo '{'
     echo '  echo api="$api"'
     echo '  echo patch="$patch"'
     echo '  echo head="$head"'
     echo '  echo'
     echo '  git ls-tree -r "$head"'
-    echo '} > "$outputDir/$manifest0"'
+    echo '} > "$outputDir/${manifestFile:?}"'
+    echo 'unset manifestFile'
 fi
 
 #------------------------------------------------------------------------------
-# Add in modules
+# OpenFOAM modules/plugins
 
 # Recursive addition of submodule content.
 # NB: must be called from within the respective parent directory.
@@ -411,11 +474,15 @@ fi
 #
 #    packModule abc  (implied cd)
 #       packModule abc/def
-#          packModule abc/def/hij
+#          packModule abc/def/ghi
 #
 packModule()
 {
     local parent="$1"
+    local filter="$2"
+
+    # Using filter=true means accept everything
+    if [ "$filter" = true ]; then unset filter; fi
 
     (
         if [ -n "$parent" ]
@@ -424,34 +491,51 @@ packModule()
         fi
 
         git ls-tree -d HEAD | \
-        while read mode gittype sha1 module
+        while read mode gittype sha1 moduleName
         do
             [ "$gittype" == commit ] || continue
 
-            case "$module" in
+            case "$moduleName" in
             (. | ./)
                 echo
                 echo "# ----"
-                echo "# module $parent : not initialized?"
+                echo "# submodule $parent : not initialized?"
                 echo "# ----"
                 continue
                 ;;
             esac
 
             # Fully qualified
-            module="$parent${parent:+/}$module"
+            module="$parent${parent:+/}$moduleName"
+            moduleName="${moduleName##*/}"
 
             echo
-            echo "# module"
+            echo "# submodule"
             echo "module='$module'"
             echo "commit='$sha1'"
-            echo "tmpTarFile='$tarName-${module##*/}.tar'"
+
+            # Simplistic filtering
+            if [ -n "$filter" ] && [ "${filter/$moduleName/}" = "$filter" ]
+            then
+                echo "# ----"
+                echo '{'
+                echo '  echo'
+                echo '  echo "$module"'
+                echo '  echo commit="$commit"'
+                echo '  echo "# not exported"'
+                echo '  echo'
+                echo '} >> "$outputDir/${manifestFile:?}"'
+                continue
+            fi
+
+            # Intermediate tar file for module contents
+            echo "tmpTarFile='$tarName-$moduleName.tar'"
             echo "# ----"
             echo '('
             echo '  cd "$module" || exit'
             echo '  newPrefix="$prefixDir${prefixDir:+/}$module"'
             echo '  git -c tar.umask=0022 archive --format=tar --prefix="$newPrefix/" -o "$outputDir/$tmpTarFile" "$commit" || exit'
-            # Without {test,tests,validation} directories (potentially large)
+            echo '  # Without {test,tests,validation} directories (potentially large)'
             echo '  tar --delete -f "$outputDir/$tmpTarFile" "$newPrefix/test" "$newPrefix/tests" "$newPrefix/validation" 2>/dev/null'
             echo '  tar -Af "$outputDir/$tarName.tar" "$outputDir/$tmpTarFile"'
             echo '  {'
@@ -461,32 +545,58 @@ packModule()
             echo '    echo'
             echo '    # Without {test,tests,validation} directories'
             echo '    git ls-tree -r "$commit" | sed -e '"'"'/\t\(test\|\tests\|validation\)\//d'"'"
-            echo '  } >> "$outputDir/$manifest1"'
+            echo '  } >> "$outputDir/${manifestFile:?}"'
             echo ')'
-            # Remove tmp tar file, even if git archive failed (creates zero-sized tar)
+            echo "# ----"
             echo 'rm -f "$outputDir/$tmpTarFile"'
 
+            # No filter for lower levels...
             packModule "$module"
         done
     )
 }
 
 
-if [ "$withModules" != false ]
+# modules/
+if [ "$select_modules" != false ]
 then
     echo
-    echo '# Modules'
+    echo '# modules/'
+    echo 'manifestFile="$manifest1"'
     echo '{'
     echo '  echo "# OpenFOAM modules"'
     echo '  echo api="$api"'
     echo '  echo patch="$patch"'
     echo '  echo head="$head"'
-    echo '} > "$outputDir/$manifest1"'
+    echo '} > "$outputDir/${manifestFile:?}"'
 
-    packModule modules
+    # With all or specified modules
+    packModule modules "$select_modules"
 
     echo
-    echo '{ echo; echo "# End"; } >> "$outputDir/$manifest1"'
+    echo '{ echo; echo "# End"; } >> "$outputDir/${manifestFile:?}"'
+    echo 'unset manifestFile'
+fi
+
+# plugins/
+if [ "$select_plugins" != false ]
+then
+    echo
+    echo '# plugins/'
+    echo 'manifestFile="$manifest2"'
+    echo '{'
+    echo '  echo "# OpenFOAM plugins"'
+    echo '  echo api="$api"'
+    echo '  echo patch="$patch"'
+    echo '  echo head="$head"'
+    echo '} > "$outputDir/${manifestFile:?}"'
+
+    # With all or specified plugins
+    packModule plugins "$select_plugins"
+
+    echo
+    echo '{ echo; echo "# End"; } >> "$outputDir/${manifestFile:?}"'
+    echo 'unset manifestFile'
 fi
 
 #------------------------------------------------------------------------------
@@ -494,10 +604,10 @@ fi
 # Decode '@' in the names as '/' directory separator
 
 echo
-echo "echo 'Adding build-info and manifest files'"
+echo echo 'Adding build-info and manifest files'
 echo 'if pushd "$outputDir"; then'
-echo "tar --owner=root --group=root --append --transform='s|@|/|g' -v -f \"\$tarName.tar\" \"\$buildInfo\" \"\$manifest0\" \"\$manifest1\""
-echo 'rm -f "$buildInfo" "$manifest0" "$manifest1"'
+echo "tar --owner=root --group=root --append --transform='s|@|/|g' -v -f \"\$tarName.tar\" \"\$buildInfo\" \"\$manifest0\" \"\$manifest1\" \"\$manifest2\""
+echo 'rm -f "$buildInfo" "$manifest0" "$manifest1" "$manifest2"'
 echo 'popd; fi'
 
 echo
@@ -509,18 +619,18 @@ echo "# -----------------------"
 # Compression
 
 case "$compress" in
-    ('')
+('')
     echo "No compression requested" 1>&2
     ;;
 
-    (gz | gzip)
+(gz | gzip)
     echo "Using gzip compression" 1>&2
     echo 'gzip -f -9 "$outputDir/$tarName.tar"'
     echo
     echo '# End of compression'
     ;;
 
-    (tgz)
+(tgz)
     echo "Using gzip compression with tgz ending" 1>&2
     echo 'gzip -c -9 "$outputDir/$tarName.tar" > "$outputDir/$tarName.tgz"'
     echo 'rm -f "$outputDir/$tarName.tar"'
@@ -528,28 +638,28 @@ case "$compress" in
     echo '# End of compression'
     ;;
 
-    (bz | bzip | bzip2)
+(bz | bzip | bzip2)
     echo "Using bzip2 compression" 1>&2
     echo 'bzip2 -f -9 "$outputDir/$tarName.tar"'
     echo
     echo '# End of compression'
     ;;
 
-    (xz)
+(xz)
     echo "Using xz compression" 1>&2
     echo 'xz -f -9 "$outputDir/$tarName.tar"'
     echo
     echo '# End of compression'
     ;;
 
-    (zst | zstd)
+(zst | zstd)
     echo "Using zstd compression" 1>&2
     echo 'zstd --rm -f -9 "$outputDir/$tarName.tar"'
     echo
     echo '# End of compression'
     ;;
 
-    (*)
+(*)
     echo "Unknown compression scheme: $compress" 1>&2
     ;;
 esac

doc/Build.md

@@ -1,23 +1,23 @@
 ## Getting the code
 
-Links to all code packs are available on https://dl.openfoam.com. For OpenFOAM-v2312:
+Links to all code packs are available on https://dl.openfoam.com. For OpenFOAM-v2406:
 
 - https://dl.openfoam.com/source/latest/
-- Source: https://dl.openfoam.com/source/v2312/OpenFOAM-v2312.tgz
-- ThirdParty: https://dl.openfoam.com/source/v2312/ThirdParty-v2312.tgz
+- Source: https://dl.openfoam.com/source/v2406/OpenFOAM-v2406.tgz
+- ThirdParty: https://dl.openfoam.com/source/v2406/ThirdParty-v2406.tgz
 
 ## OpenFOAM® Quick Build Guide
 
 Prior to building, ensure that the [system requirements][link openfoam-require]
 are satisfied (including any special [cross-compiling][wiki-cross-compile]
 considerations), and source the correct OpenFOAM environment.
-For example, for the OpenFOAM-v2312 version:
+For example, for the OpenFOAM-v2406 version:
 ```
-source <installation path>/OpenFOAM-v2312/etc/bashrc
+source <installation path>/OpenFOAM-v2406/etc/bashrc
 ```
 e.g. if installed under the `~/openfoam` directory
 ```
-source ~/openfoam/OpenFOAM-v2312/etc/bashrc
+source ~/openfoam/OpenFOAM-v2406/etc/bashrc
 ```
 
 
@@ -157,4 +157,4 @@ More details in the [ThirdParty build guide][link third-build].
 [wiki-config]: https://develop.openfoam.com/Development/openfoam/-/wikis/configuring
 
 ---
-Copyright 2019-2023 OpenCFD Ltd
+Copyright 2019-2024 OpenCFD Ltd

doc/BuildIssues.md

@@ -1,4 +1,4 @@
-## Known Build Issues (v2312, v2306, v2212, v2206, v2112, v2106, v2012)
+## Known Build Issues (v2406, v2312, v2306, v2212, v2206, v2112, v2106, v2012)
 
 ### Windows cross-compilation
 
@@ -188,4 +188,4 @@ and attempt to install a `paraview~qt` version instead.
 [wiki-config]: https://develop.openfoam.com/Development/openfoam/-/wikis/configuring
 
 ---
-Copyright 2019-2023 OpenCFD Ltd
+Copyright 2019-2024 OpenCFD Ltd

etc/bashrc

@@ -6,7 +6,7 @@
 #    \\/     M anipulation  |
 #------------------------------------------------------------------------------
 #     Copyright (C) 2011-2016 OpenFOAM Foundation
-#     Copyright (C) 2016-2023 OpenCFD Ltd.
+#     Copyright (C) 2016-2024 OpenCFD Ltd.
 #------------------------------------------------------------------------------
 # License
 #     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -55,7 +55,7 @@
 
 # [WM_PROJECT_VERSION] - A human-readable version name
 # A development version is often named 'com' - as in www.openfoam.com
-export WM_PROJECT_VERSION=v2312
+export WM_PROJECT_VERSION=v2406
 
 #------------------------------------------------------------------------------
 # Configuration environment variables.

etc/caseDicts/annotated/topoSetSourcesDict

@@ -34,6 +34,16 @@ cellSet_doc
         // or
         min     (0 0 0);
         max     (1 1 1);
+
+        // Optional the cell centres can be transformed before being used. This
+        // can be used to simulate moving selections in transient simulations.
+        // Any solidBodyMotionFunction can be used. Note: the same optional
+        // transformation can be used for any selection requiring face centres
+        // or point locations.
+        solidBodyMotionFunction linearMotion
+        {
+            velocity    (-0.1 -0.1 0);
+        }
     }
 
 
@@ -58,6 +68,8 @@ cellSet_doc
         radius  5.0;
         // optional
         innerRadius 1.0;
+
+        // See boxToCell for optional transformation of the cell centres
     }
 
     //- Select based on faceSet
@@ -117,6 +129,8 @@ cellSet_doc
     {
         source  nearestToCell;
         points  ((0 0 0) (1 1 1)(2 2 2));
+
+        // See boxToCell for optional transformation of the cell centres
     }
 
 
@@ -162,6 +176,8 @@ cellSet_doc
         i       (0.2 0.2 0);
         j       (-0.2 0.2 0);
         k       (10 10 10);
+
+        // See boxToCell for optional transformation of the cell centres
     }
 
 
@@ -170,6 +186,8 @@ cellSet_doc
         source       searchableSurfaceToCell;
         surfaceType  triSurfaceMesh;  // Type of surface
         surfaceName  blob.obj;        // Name for surface object
+
+        // See boxToCell for optional transformation of the cell centres
     }
 
 
@@ -188,6 +206,8 @@ cellSet_doc
 
         // optional
         innerRadius 1.0;
+
+        // See boxToCell for optional transformation of the cell centres
     }
 
 
@@ -209,6 +229,7 @@ cellSet_doc
         useSurfaceOrientation false;  // use closed surface inside/outside
                                       // test (ignores includeCut,
                                       // outsidePoints)
+        // See boxToCell for optional transformation of the cell centres
     }
 
 
@@ -257,6 +278,8 @@ faceSet_doc
         // or
         min     (0 0 0);
         max     (1 1 1);
+
+        // See boxToCell for optional transformation of the face centres
     }
 
 
@@ -286,6 +309,8 @@ faceSet_doc
         radius  5.0;
         // optional
         innerRadius 1.0;
+
+        // See boxToCell for optional transformation of the face centres
     }
 
 
@@ -352,6 +377,8 @@ faceSet_doc
         source       searchableSurfaceToFace;
         surfaceType  triSurfaceMesh;  // Type of surface
         surfaceName  blob.obj;        // Name for surface object
+
+        // See boxToCell for optional transformation of the face centres
     }
 
 
@@ -362,6 +389,8 @@ faceSet_doc
         radius  5.0;
         // optional
         innerRadius 1.0;
+
+        // See boxToCell for optional transformation of the face centres
     }
 
 
@@ -410,6 +439,8 @@ pointSet_doc
         // or
         min     (0 0 0);
         max     (1 1 1);
+
+        // See boxToCell for optional transformation of the points
     }
 
 
@@ -436,6 +467,7 @@ pointSet_doc
         radius  5.0;
         // optional
         innerRadius 1.0;
+        // See boxToCell for optional transformation of the points
     }
 
 
@@ -465,6 +497,7 @@ pointSet_doc
     {
         source  nearestToPoint;
         points  ((0 0 0) (1 1 1));
+        // See boxToCell for optional transformation of the points
     }
 
 
@@ -486,6 +519,7 @@ pointSet_doc
         source       searchableSurfaceToCell;
         surfaceType  triSurfaceMesh;  // Type of surface
         surfaceName  blob.obj;        // Name for surface object
+        // See boxToCell for optional transformation of the points
     }
 
 
@@ -496,6 +530,7 @@ pointSet_doc
         radius  5.0;
         // optional
         innerRadius 1.0;
+        // See boxToCell for optional transformation of the points
     }
 
 
@@ -508,6 +543,7 @@ pointSet_doc
                                 // (requires closed surface with consistent
                                 //  normals)
         includeOutside  false;  //   ,,    outside  ,,
+        // See boxToCell for optional transformation of the points
     }
 
 

etc/caseDicts/setConstraintTypes

@@ -14,6 +14,7 @@ cyclic
 cyclicAMI
 {
     type  cyclicAMI;
+    value $internalField;
 }
 
 cyclicACMI

etc/config.csh/adios2

@@ -5,7 +5,7 @@
 #   \\  /    A nd           | www.openfoam.com
 #    \\/     M anipulation  |
 #------------------------------------------------------------------------------
-#     Copyright (C) 2017-2023 OpenCFD Ltd.
+#     Copyright (C) 2017-2024 OpenCFD Ltd.
 #------------------------------------------------------------------------------
 # License
 #     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -20,7 +20,7 @@
 #------------------------------------------------------------------------------
 # USER EDITABLE PART: Changes made here may be lost with the next upgrade
 
-set adios2_version=ADIOS2-2.9.2
+set adios2_version=ADIOS2-2.10.1
 setenv ADIOS2_ARCH_PATH "$WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER/$adios2_version"
 
 # END OF (NORMAL) USER EDITABLE PART

etc/config.csh/hypre

@@ -5,7 +5,7 @@
 #   \\  /    A nd           | www.openfoam.com
 #    \\/     M anipulation  |
 #------------------------------------------------------------------------------
-#     Copyright (C) 2018-2021 OpenCFD Ltd.
+#     Copyright (C) 2018-2024 OpenCFD Ltd.
 #------------------------------------------------------------------------------
 # License
 #     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -30,7 +30,7 @@
 #------------------------------------------------------------------------------
 # USER EDITABLE PART: Changes made here may be lost with the next upgrade
 
-set hypre_version=hypre-2.23.0
+set hypre_version=hypre-2.31.0
 setenv HYPRE_ARCH_PATH "$WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER$WM_PRECISION_OPTION$WM_LABEL_OPTION/$hypre_version"
 
 # END OF (NORMAL) USER EDITABLE PART

etc/config.csh/paraview

@@ -6,7 +6,7 @@
 #    \\/     M anipulation  |
 #------------------------------------------------------------------------------
 #     Copyright (C) 2011-2016 OpenFOAM Foundation
-#     Copyright (C) 2016-2023 OpenCFD Ltd.
+#     Copyright (C) 2016-2024 OpenCFD Ltd.
 #------------------------------------------------------------------------------
 # License
 #     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -52,7 +52,7 @@
 #------------------------------------------------------------------------------
 # USER EDITABLE PART: Changes made here may be lost with the next upgrade
 
-set ParaView_VERSION=5.11.2
+set ParaView_VERSION=5.12.1
 set ParaView_QT=qt-system
 
 # END OF (NORMAL) USER EDITABLE PART

etc/config.csh/petsc

@@ -5,7 +5,7 @@
 #   \\  /    A nd           | www.openfoam.com
 #    \\/     M anipulation  |
 #------------------------------------------------------------------------------
-#     Copyright (C) 2018-2023 OpenCFD Ltd.
+#     Copyright (C) 2018-2024 OpenCFD Ltd.
 #------------------------------------------------------------------------------
 # License
 #     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -30,7 +30,7 @@
 #------------------------------------------------------------------------------
 # USER EDITABLE PART: Changes made here may be lost with the next upgrade
 
-set petsc_version=petsc-3.19.2
+set petsc_version=petsc-3.21.2
 setenv PETSC_ARCH_PATH "$WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER$WM_PRECISION_OPTION$WM_LABEL_OPTION/$petsc_version"
 
 # END OF (NORMAL) USER EDITABLE PART

etc/config.sh/adios2

@@ -5,7 +5,7 @@
 #   \\  /    A nd           | www.openfoam.com
 #    \\/     M anipulation  |
 #------------------------------------------------------------------------------
-#     Copyright (C) 2017-2023 OpenCFD Ltd.
+#     Copyright (C) 2017-2024 OpenCFD Ltd.
 #------------------------------------------------------------------------------
 # License
 #     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -21,7 +21,7 @@
 #------------------------------------------------------------------------------
 # USER EDITABLE PART: Changes made here may be lost with the next upgrade
 
-adios2_version=ADIOS2-2.9.2
+adios2_version=ADIOS2-2.10.1
 export ADIOS2_ARCH_PATH="$WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER/$adios2_version"
 
 # END OF (NORMAL) USER EDITABLE PART

etc/config.sh/hdf5

@@ -0,0 +1,61 @@
+#----------------------------------*-sh-*--------------------------------------
+# =========                 |
+# \\      /  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.
+#
+# File
+#     etc/config.sh/hdf5
+#     [optional/experimental: not sourced by OpenFOAM-*/etc/bashrc]
+#
+# Description
+#     Setup for HDF5 include/libraries (usually ThirdParty installation).
+#
+#     To disable its use:               hdf5_version=hdf5-none
+#     For system-wide installations:    hdf5_version=hdf5-system
+#
+#     For central installations not located under ThirdParty:
+#       1. use hdf5-system
+#       2. and provide full path for HDF5_ARCH_PATH
+#
+#------------------------------------------------------------------------------
+# USER EDITABLE PART: Changes made here may be lost with the next upgrade
+
+hdf5_version=hdf5-1.14.4.3
+export HDF5_ARCH_PATH="$WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER/$hdf5_version"
+
+# END OF (NORMAL) USER EDITABLE PART
+#------------------------------------------------------------------------------
+if [ -n "$FOAM_VERBOSE" ] && [ -n "$PS1" ]
+then
+    echo "Using hdf5 ($hdf5_version)  ->  $HDF5_ARCH_PATH" 1>&2
+fi
+
+if command -v _foamAddLibAuto >/dev/null
+then
+    # Normal sourcing (not makeHDF5)
+
+    _foamAddLibAuto $HDF5_ARCH_PATH
+
+    unset hdf5_version
+
+elif [ "$1" = "-force" ]
+then
+    # Forced command-line sourcing
+
+    if output="$($WM_PROJECT_DIR/bin/tools/lib-dir -sh $HDF5_ARCH_PATH 2>/dev/null)"
+    then
+        eval "$output"
+    fi
+
+    unset hdf5_version output
+fi
+
+
+#------------------------------------------------------------------------------

etc/config.sh/hypre

@@ -5,7 +5,7 @@
 #   \\  /    A nd           | www.openfoam.com
 #    \\/     M anipulation  |
 #------------------------------------------------------------------------------
-#     Copyright (C) 2018-2021 OpenCFD Ltd.
+#     Copyright (C) 2018-2024 OpenCFD Ltd.
 #------------------------------------------------------------------------------
 # License
 #     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -27,7 +27,7 @@
 #------------------------------------------------------------------------------
 # USER EDITABLE PART: Changes made here may be lost with the next upgrade
 
-hypre_version=hypre-2.23.0
+hypre_version=hypre-2.31.0
 export HYPRE_ARCH_PATH="$WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER$WM_PRECISION_OPTION$WM_LABEL_OPTION/$hypre_version"
 
 # END OF (NORMAL) USER EDITABLE PART

etc/config.sh/paraview

@@ -6,7 +6,7 @@
 #    \\/     M anipulation  |
 #------------------------------------------------------------------------------
 #     Copyright (C) 2011-2016 OpenFOAM Foundation
-#     Copyright (C) 2016-2023 OpenCFD Ltd.
+#     Copyright (C) 2016-2024 OpenCFD Ltd.
 #------------------------------------------------------------------------------
 # License
 #     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -57,7 +57,7 @@
 #------------------------------------------------------------------------------
 # USER EDITABLE PART: Changes made here may be lost with the next upgrade
 
-ParaView_VERSION=5.11.2
+ParaView_VERSION=5.12.1
 ParaView_QT=qt-system
 
 # END OF (NORMAL) USER EDITABLE PART

etc/config.sh/petsc

@@ -5,7 +5,7 @@
 #   \\  /    A nd           | www.openfoam.com
 #    \\/     M anipulation  |
 #------------------------------------------------------------------------------
-#     Copyright (C) 2018-2023 OpenCFD Ltd.
+#     Copyright (C) 2018-2024 OpenCFD Ltd.
 #------------------------------------------------------------------------------
 # License
 #     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -27,7 +27,7 @@
 #------------------------------------------------------------------------------
 # USER EDITABLE PART: Changes made here may be lost with the next upgrade
 
-petsc_version=petsc-3.19.2
+petsc_version=petsc-3.21.2
 export PETSC_ARCH_PATH="$WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER$WM_PRECISION_OPTION$WM_LABEL_OPTION/$petsc_version"
 
 # END OF (NORMAL) USER EDITABLE PART

etc/cshrc

@@ -6,7 +6,7 @@
 #    \\/     M anipulation  |
 #------------------------------------------------------------------------------
 #     Copyright (C) 2011-2016 OpenFOAM Foundation
-#     Copyright (C) 2016-2023 OpenCFD Ltd.
+#     Copyright (C) 2016-2024 OpenCFD Ltd.
 #------------------------------------------------------------------------------
 # License
 #     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -55,7 +55,7 @@
 
 # [WM_PROJECT_VERSION] - A human-readable version name
 # A development version is often named 'com' - as in www.openfoam.com
-setenv WM_PROJECT_VERSION v2312
+setenv WM_PROJECT_VERSION v2406
 
 #------------------------------------------------------------------------------
 # Configuration environment variables.

modules/README.md

@@ -1,12 +1,17 @@
 [[_TOC_]]
 
-## OpenFOAM Modules
+## OpenFOAM modules
 
 This directory is a location for additional OpenFOAM components or
 tools to placed and have them built as part of the normal OpenFOAM
 build process. It is assumed that each subdirectory contain an
 appropriate `Allwmake` (or `Allwmake.override`) file.
 
+The primary distinction between `modules` and `plugins` is that `modules` are
+mainly maintained and released by OpenFOAM developers, whereas `plugins` are an
+open and welcoming area for add-ons, predominantly maintained and driven by
+OpenFOAM community members and groups.
+
 ### How to use
 
 On the first use, you need to register the submodules, and then update them.
@@ -91,7 +96,7 @@ To build the doxygen information for the components, it is also
 necessary to link the directories to the doc/ subdirectory.
 This is a purely manual operation.
 
-### Developer Information
+### Developer information
 
 #### Build locations
 

plugins/Allwmake

@@ -0,0 +1,42 @@
+#!/bin/sh
+cd "${0%/*}" || exit                            # Run from this directory
+targetType=libso
+. "${WM_PROJECT_DIR:?}"/wmake/scripts/AllwmakeParseArguments
+
+#------------------------------------------------------------------------------
+# Default build into OpenFOAM project locations unless specified with
+# -prefix or FOAM_MODULE_PREFIX env varable
+
+# Long form to avoid dash 0.5.8 error (issue #1757)
+[ -n "$FOAM_MODULE_PREFIX" ] || FOAM_MODULE_PREFIX="${FOAM_LIBBIN%/*}"
+export FOAM_MODULE_PREFIX
+
+echo "========================================"
+case "$FOAM_MODULE_PREFIX" in
+(false | none)
+    echo "OpenFOAM modules disabled (prefix=${FOAM_MODULE_PREFIX})"
+    echo
+    exit 0
+esac
+
+echo "prefix = $FOAM_MODULE_PREFIX"
+echo
+echo "    ignoring possible compilation errors"
+echo "    make certain to check the output file"
+echo
+set +e
+export WM_CONTINUE_ON_ERROR=true
+
+#------------------------------------------------------------------------------
+
+for name in $(./list-plugins)
+do
+    if [ -d "$name" ]
+    then
+        # Use wmake -all instead of Allwmake to allow for overrides
+        ( cd "$name" && wmake -all $targetType )
+    fi
+done
+
+
+#------------------------------------------------------------------------------

plugins/README.md

@@ -0,0 +1,139 @@
+[[_TOC_]]
+
+## OpenFOAM plugins
+
+This directory is a location for additional OpenFOAM components or
+tools to placed and have them built as part of the normal OpenFOAM
+build process. It is assumed that each subdirectory contain an
+appropriate `Allwmake` (or `Allwmake.override`) file.
+
+The primary distinction between `modules` and `plugins` is that `modules` are
+mainly maintained and released by OpenFOAM developers, whereas `plugins` are an
+open and welcoming area for add-ons, predominantly maintained and driven by
+OpenFOAM community members and groups.
+
+### How to use
+
+On the first use, you need to register the submodules, and then update them.
+You can execute both steps for all the available submodules (including the
+nested ones) as follows while you are at `$WM_PROJECT_DIR`:
+
+```bash
+cd $WM_PROJECT_DIR
+
+git submodule update --init --recursive
+```
+
+Executing this single-line command clones all the submodules from their
+respective repositories and prepares them for compilation. Note that you can
+also make only a certain group of submodules ready by explicitly specifying the
+requested submodules' names at the end of the command above. For example, if
+you would like to use only the `turbulence-community` submodule, you specify:
+
+```bash
+git submodule update --init --recursive plugins/turbulence-community
+```
+
+You can display information about the status of submodules as follows:
+
+```bash
+git submodule status --recursive
+```
+
+An easy way to see which submodules are actually in use:
+
+```bash
+cat .gitmodules
+```
+
+Which will reveal content resembling the following:
+```
+[submodule "xyz"]
+    path = plugins/xyz
+    url = ...
+...
+```
+
+If you need to remove a specific submodule or wish to restart the process,
+you can simply carry out the task as follows:
+
+```bash
+git submodule deinit plugins/turbulence-community
+```
+
+This command deregisters the specified submodule and clears the
+`plugins/turbulence-community` directory.
+
+A quick overview of `git submodules` can be found in this
+[*blog*][blog git-submodule] with full details in the
+[*manpage*][man git-submodule].
+
+### Build locations
+
+Any individual _plugin_ will normally also be able to exist outside of
+the plugins directory structure and will typically build into user
+locations (`$FOAM_USER_APPBIN`, `$FOAM_USER_LIBBIN`).
+
+When compiled from the top-level OpenFOAM `Allwmake` or the
+`plugins/Allwmake`, they should build into OpenFOAM project locations
+(`$FOAM_APPBIN`, `$FOAM_LIBBIN`). This can be adjusted by
+supplying an alternative `-prefix=` to the corresponding Allwmake
+command.
+
+| Command    | Install location |
+|------------|------------------|
+| ./Allwmake -prefix=user | `$FOAM_USER_APPBIN`, `$FOAM_USER_LIBBIN` |
+| ./Allwmake -prefix=group | `$FOAM_SITE_APPBIN`, `$FOAM_SITE_LIBBIN` |
+| ./Allwmake -prefix=openfoam | `$FOAM_APPBIN`, `$FOAM_LIBBIN` |
+| ./Allwmake -prefix=/some/pathname | `/some/pathname/bin`, `/some/pathname/lib` |
+
+### Documentation (doxygen)
+
+To build the doxygen information for the components, it is also
+necessary to link the directories to the doc/ subdirectory.
+This is a purely manual operation.
+
+### Developer information
+
+#### Build locations
+
+To accomodate building into various locations, the plugin code should
+be adapted with the following changes:
+
+- ***Make/files***
+   ```
+   ...
+   EXE = $(FOAM_MODULE_APPBIN)/someExecutable
+
+   LIB = $(FOAM_MODULE_LIBBIN)/libSomeLibrary
+   ```
+
+- `Make/options` should include this
+  ```
+  include $(GENERAL_RULES)/module-path-user
+  ...
+  ```
+
+The following changes to `Make/options` are universally applicable
+(ie, work with older or other versions of OpenFOAM), but more verbose.
+
+- `Make/options` with the following
+  ```
+  sinclude $(GENERAL_RULES)/module-path-user
+
+  /* Failsafe - user locations */
+  ifeq (,$(FOAM_MODULE_APPBIN))
+  FOAM_MODULE_APPBIN = $(FOAM_USER_APPBIN)
+  endif
+  ifeq (,$(FOAM_MODULE_LIBBIN))
+  FOAM_MODULE_LIBBIN = $(FOAM_USER_LIBBIN)
+  endif
+  ...
+  ```
+
+<!-- General Information -->
+
+[man git-submodule]:  https://git-scm.com/docs/git-submodule
+[blog git-submodule]: http://blog.joncairns.com/2011/10/how-to-use-git-submodules/
+
+---

plugins/list-plugins

@@ -0,0 +1,85 @@
+#!/bin/sh
+#------------------------------------------------------------------------------
+# =========                 |
+# \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+#  \\    /   O peration     |
+#   \\  /    A nd           | www.openfoam.com
+#    \\/     M anipulation  |
+#------------------------------------------------------------------------------
+#     Copyright (C) 2020 OpenCFD Ltd.
+#------------------------------------------------------------------------------
+# SPDX-License-Identifier: (GPL-3.0-or-later)
+#
+# Script
+#     list-modules
+#
+# Description
+#     List module directories
+#     - each first-level directory with an Allwmake file
+#
+#------------------------------------------------------------------------------
+cd "${0%/*}" || exit                            # Run from this directory
+
+printHelp() {
+    cat<< HELP 1>&2
+
+Usage: ${0##*/} [OPTION]
+options:
+  -help             Display help and exit
+
+List module directories - each first-level directory with an Allwmake file
+
+HELP
+
+    exit 0  # A clean exit
+}
+
+# Report error and exit
+die()
+{
+    exec 1>&2
+    echo
+    echo "Error encountered:"
+    while [ "$#" -ge 1 ]; do echo "    $1"; shift; done
+    echo
+    echo "See '${0##*/} -help' for usage"
+    echo
+    exit 1
+}
+
+#------------------------------------------------------------------------------
+
+# Parse options
+while [ "$#" -gt 0 ]
+do
+    case "$1" in
+    -h | -help*)  # Short help
+        printHelp
+        ;;
+
+    *)
+        die "Unknown option/argument: '$1'"
+        ;;
+    esac
+    shift
+done
+
+
+# Each first-level directory with an Allwmake file
+for moduleName in *
+do
+    if [ -f "$moduleName/Allwmake" ]
+    then
+        case "$moduleName" in
+        # Skip some directory names
+        (build | doc | platform*)
+            ;;
+        (*)
+            echo "$moduleName"
+            ;;
+        esac
+    fi
+done
+
+
+#------------------------------------------------------------------------------

src/OpenFOAM/Make/files

@@ -185,6 +185,7 @@ $(expr)/scanToken/exprScanToken.C
 
 $(expr)/traits/exprTraits.C
 $(expr)/value/exprValue.C
+$(expr)/value/exprValueFieldTag.C
 
 $(expr)/exprDriver/exprDriver.C
 $(expr)/exprDriver/exprDriverFields.C
@@ -445,7 +446,6 @@ $(lduMatrix)/solvers/PBiCGStab/PBiCGStab.C
 $(lduMatrix)/solvers/FPCG/FPCG.C
 $(lduMatrix)/solvers/PPCG/PPCG.C
 $(lduMatrix)/solvers/PPCR/PPCR.C
-$(lduMatrix)/solvers/GMRES/gmresSolver.C
 
 $(lduMatrix)/smoothers/GaussSeidel/GaussSeidelSmoother.C
 $(lduMatrix)/smoothers/symGaussSeidel/symGaussSeidelSmoother.C

src/OpenFOAM/expressions/exprResult/exprResult.H

@@ -146,11 +146,11 @@ class exprResult
         //- Type-checked determination of centre value (min/max)
         //  \return True if the type check was satisfied
         template<class Type>
-        bool setAverageValueChecked(const bool parRun = Pstream::parRun());
+        bool setAverageValueChecked(const bool parRun = UPstream::parRun());
 
         //- Type-checked determination of average bool value
         //  \return True if the type check was satisfied
-        bool setAverageValueCheckedBool(const bool parRun = Pstream::parRun());
+        bool setAverageValueCheckedBool(const bool parRun = UPstream::parRun());
 
         //- Type-checked copy of field
         //  \return True if the type check was satisfied
@@ -385,7 +385,7 @@ public:
         //- Test if field corresponds to a single-value and thus uniform.
         //  Uses field min/max to establish uniformity.
         //  Test afterwards with isUniform()
-        void testIfSingleValue(const bool parRun = Pstream::parRun());
+        void testIfSingleValue(const bool parRun = UPstream::parRun());
 
 
     // Set results
@@ -437,7 +437,7 @@ public:
         (
             const label size,
             const bool noWarn,
-            const bool parRun = Pstream::parRun()
+            const bool parRun = UPstream::parRun()
         ) const;
 
         //- Get a reduced result

src/OpenFOAM/expressions/value/exprValue.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021-2023 OpenCFD Ltd.
+    Copyright (C) 2021-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -26,24 +26,42 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "exprValue.H"
+#include "error.H"
 #include "ITstream.H"
 #include "Switch.H"
 #include <cstring>  // For memcpy, memset
 
+// * * * * * * * * * * * * * * * *  Details  * * * * * * * * * * * * * * * * //
+
+void Foam::expressions::Detail::exprValueUnion::notSpecialized
+(
+    const std::string& msg
+) noexcept
+{
+    FatalErrorInFunction
+        << "non-specialized: " << msg.c_str() << endl
+        << abort(FatalError);
+}
+
+
 // * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
+#if 0
+// General way to add tokens for VectorSpace types
+// (caller excludes none/invalid)
 template<class Type>
-static void fillTokens(const Type& val, tokenList& toks)
+static void addTokens(tokenList& toks, const Type& val)
 {
     const direction nCmpt = pTraits<Type>::nComponents;
     const direction nParen = 2*(pTraits<Type>::rank || (nCmpt > 1) ? 1 : 0);
 
-    toks.resize_nocopy(nCmpt + nParen);
+    const label nOld = toks.size();
+    toks.resize(nOld + label(nCmpt + nParen));
 
-    auto iter = toks.begin();
+    auto iter = toks.begin(nOld);
 
     if (nParen)
     {
@@ -67,28 +85,57 @@ static void fillTokens(const Type& val, tokenList& toks)
 
 //- Specialized for bool
 template<>
-void fillTokens<bool>(const bool& val, tokenList& toks)
+void addTokens<bool>(tokenList& toks, const bool& val)
 {
-    toks.resize_nocopy(1);
-    toks.front() = token::boolean(val);
+    toks.emplace_back() = token::boolean(val);
 }
 
 
 //- Specialized for label
 template<>
-void fillTokens<label>(const label& val, tokenList& toks)
+void addTokens<label>(tokenList& toks, const label& val)
 {
-    toks.resize_nocopy(1);
-    toks.front() = val;
+    toks.emplace_back() = val;
 }
 
 
 //- Specialized for scalar
 template<>
-void fillTokens<scalar>(const scalar& val, tokenList& toks)
+void addTokens<scalar>(tokenList& toks, const scalar& val)
 {
-    toks.resize_nocopy(1);
-    toks.front() = val;
+    toks.emplace_back() = val;
+}
+
+#endif
+
+
+// General output of type (caller excludes none/invalid)
+template<class Type>
+static void putType(Ostream& os, const Type& val)
+{
+    os  << val;
+}
+
+
+//- Specialized for bool.
+//- Write as (true/false) via Switch to avoid bool/label ambiguity
+template<>
+void putType<bool>(Ostream& os, const bool& val)
+{
+    // Note: prefer Switch() vs (std::ios::boolalpha) to avoid any
+    // potential locale issues.
+    os  << Switch(val);
+}
+
+
+//- Specialized for scalar.
+//- Write with '.' to avoid scalar/label ambiguity
+template<>
+void putType<scalar>(Ostream& os, const scalar& val)
+{
+    const auto oldflags = os.setf(std::ios::showpoint);
+    os  << val;
+    os.flags(oldflags);  // Restore
 }
 
 } // End namespace Foam
@@ -226,6 +273,7 @@ void Foam::expressions::exprValue::deepCopy(const exprValue& rhs)
 }
 
 
+#if 0
 Foam::tokenList Foam::expressions::exprValue::tokens(bool prune) const
 {
     // Handling for NONE, INVALID:
@@ -237,24 +285,28 @@ Foam::tokenList Foam::expressions::exprValue::tokens(bool prune) const
 
     tokenList toks;
 
-    if (!prune)
-    {
-        if (typeCode_ == expressions::valueTypeCode::NONE)
-        {
-            toks.resize(2);
-            toks.front() = token::BEGIN_LIST;
-            toks.back() = token::END_LIST;
-            return toks;
-        }
-        else if (typeCode_ == expressions::valueTypeCode::INVALID)
-        {
-            toks.emplace_back(word("bad"));
-            return toks;
-        }
-    }
-
     switch (typeCode_)
     {
+        case expressions::valueTypeCode::NONE :
+        {
+            if (!prune)
+            {
+                toks.resize(2);
+                toks.front() = token::BEGIN_LIST;
+                toks.back() = token::END_LIST;
+            }
+            break;
+        }
+
+        case expressions::valueTypeCode::INVALID :
+        {
+            if (!prune)
+            {
+                toks.emplace_back(word("bad"));
+            }
+            break;
+        }
+
         #undef  doLocalCode
         #define doLocalCode(Type, UnusedParam)                        \
                                                                       \
@@ -263,7 +315,7 @@ Foam::tokenList Foam::expressions::exprValue::tokens(bool prune) const
             const Type* dataPtr = data_.get<Type>();                  \
             if (dataPtr)                                              \
             {                                                         \
-                fillTokens<Type>(*dataPtr, toks);                     \
+                addTokens<Type>(toks, *dataPtr);                      \
             }                                                         \
             break;                                                    \
         }
@@ -277,6 +329,7 @@ Foam::tokenList Foam::expressions::exprValue::tokens(bool prune) const
 
     return toks;
 }
+#endif
 
 
 void Foam::expressions::exprValue::write(Ostream& os, bool prune) const
@@ -288,22 +341,26 @@ void Foam::expressions::exprValue::write(Ostream& os, bool prune) const
     // With prune:
     //   - no output for either
 
-    if (!prune)
-    {
-        if (typeCode_ == expressions::valueTypeCode::NONE)
-        {
-            os << token::BEGIN_LIST << token::END_LIST;
-            return;
-        }
-        else if (typeCode_ == expressions::valueTypeCode::INVALID)
-        {
-            os << word("bad");
-            return;
-        }
-    }
-
     switch (typeCode_)
     {
+        case expressions::valueTypeCode::NONE :
+        {
+            if (!prune)
+            {
+                os << token::BEGIN_LIST << token::END_LIST;
+            }
+            break;
+        }
+
+        case expressions::valueTypeCode::INVALID :
+        {
+            if (!prune)
+            {
+                os << word("bad");
+            }
+            break;
+        }
+
         #undef  doLocalCode
         #define doLocalCode(Type, UnusedParam)                        \
                                                                       \
@@ -312,7 +369,7 @@ void Foam::expressions::exprValue::write(Ostream& os, bool prune) const
             const Type* dataPtr = data_.get<Type>();                  \
             if (dataPtr)                                              \
             {                                                         \
-                os << *dataPtr;                                       \
+                putType(os, *dataPtr);                                \
             }                                                         \
             break;                                                    \
         }
@@ -381,34 +438,27 @@ bool Foam::expressions::exprValue::readTokens(ITstream& is)
 
     const valueTypeCode whichCode(exprValue::peekType(is));
 
-    if (whichCode == expressions::valueTypeCode::NONE)
-    {
-        typeCode_ = whichCode;
-        is.skip(2);  // Skip tokens: '( )'
-        return true;
-    }
-
-    // This one should be rare or even impossible
-    if (whichCode == expressions::valueTypeCode::INVALID)
-    {
-        typeCode_ = whichCode;
-
-        if (is.bad())
-        {
-            return false;
-        }
-
-        const token& tok0 = is.peek();
-
-        if (tok0.isWord("bad"))
-        {
-            is.skip(1);  // Skip token: "bad"
-            return true;
-        }
-    }
-
     switch (whichCode)
     {
+        case expressions::valueTypeCode::NONE :
+        {
+            typeCode_ = whichCode;
+            is.skip(2);  // Skip tokens: '( )'
+            return true;
+        }
+
+        // This one should be rare or even impossible
+        case expressions::valueTypeCode::INVALID :
+        {
+            typeCode_ = whichCode;
+            if (!is.bad() && is.peek().isWord("bad"))
+            {
+                is.skip(1);  // Skip token: "bad"
+                return true;
+            }
+            return false;  // Some type of failure..
+        }
+
         #undef  doLocalCode
         #define doLocalCode(Type, UnusedParam)                        \
                                                                       \
@@ -430,6 +480,58 @@ bool Foam::expressions::exprValue::readTokens(ITstream& is)
 }
 
 
+int Foam::expressions::exprValue::compare
+(
+    const exprValue& rhs
+) const
+{
+    if (typeCode_ != rhs.typeCode_)
+    {
+        // First compare by type
+        return (int(typeCode_) - int(rhs.typeCode_));
+    }
+    else if ((this == &rhs) || !good())
+    {
+        // Identical: same object or not good
+        // (ie, no further comparison possible)
+        return 0;
+    }
+
+    // Types are identical (and good)
+    // - compare by value.
+    // This is even messier than usual, since can only rely on
+    // operator< being defined
+
+    switch (typeCode_)
+    {
+        #undef  doLocalCode
+        #define doLocalCode(Type, UnusedParam)                        \
+                                                                      \
+        case expressions::valueTypeCode::type_##Type :                \
+        {                                                             \
+            const Type* a = data_.get<Type>();                        \
+            const Type* b = rhs.data_.get<Type>();                    \
+            return                                                    \
+            (                                                         \
+                (a && b)                                              \
+              ? ((*a < *b) ? -1 : (*b < *a) ? 1 : 0)                  \
+              : 0                                                     \
+            );                                                        \
+            break;                                                    \
+        }
+
+        FOR_ALL_EXPR_VALUE_TYPES(doLocalCode);
+        #undef doLocalCode
+
+        // exprValue may only be a subset of valueTypeCode types
+        default: break;
+    }
+
+    // Should not happen
+    return 0;
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 bool Foam::expressions::exprValue::operator==(const exprValue& rhs) const
@@ -470,8 +572,7 @@ bool Foam::expressions::exprValue::operator==(const exprValue& rhs) const
 
 bool Foam::expressions::exprValue::operator<(const exprValue& rhs) const
 {
-    // Not yet sortable
-    return false;
+    return (this->compare(rhs) < 0);
 }
 
 

src/OpenFOAM/expressions/value/exprValue.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021-2023 OpenCFD Ltd.
+    Copyright (C) 2021-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -42,10 +42,8 @@ SourceFiles
 #ifndef Foam_expressions_exprValue_H
 #define Foam_expressions_exprValue_H
 
-#include "exprTraits.H"
-#include "error.H"
+#include "exprTraits.H"  // For valueTypeCode and label, scalar, vector etc.
 #include "contiguous.H"
-#include "tokenList.H"
 #include "InfoProxy.H"
 #include <typeinfo>
 
@@ -102,12 +100,7 @@ union exprValueUnion
         //  This seems to be the best way to catch programming errors
         //  since static_assert does not help here.
         //  The 'noexcept' is slightly misleading (needed for the caller)
-        static void notSpecialized(const std::string& msg) noexcept
-        {
-            FatalErrorInFunction
-                << "non-specialized: " << msg.c_str() << endl
-                << abort(FatalError);
-        }
+        static void notSpecialized(const std::string& msg) noexcept;
 
         //- Return read pointer to typed union member,
         //- which is nullptr for unspecialized versions
@@ -175,13 +168,8 @@ class exprValue
         //- Copy assignment
         void deepCopy(const exprValue& rhs);
 
-
 public:
 
-    //- Runtime type information
-    ClassName("exprValue");
-
-
     // Constructors
 
         //- Default construct (zero-initialized) as 'none'
@@ -244,7 +232,7 @@ public:
         //- The value type code
         expressions::valueTypeCode typeCode() const noexcept
         {
-           return typeCode_;
+            return typeCode_;
         }
 
         //- True if the value type is not none/invalid
@@ -277,6 +265,9 @@ public:
         //  \return True on success
         bool readTokens(ITstream& is);
 
+        //- Compare (type,value)
+        int compare(const exprValue& rhs) const;
+
 
     // Typed Access
 
@@ -350,7 +341,7 @@ public:
         //- Compare (type,value) for inequality
         bool operator!=(const exprValue& rhs) const { return !(*this == rhs); }
 
-        //- Compare (type,value) - currently not implemented.
+        //- Compare (type,value)
         bool operator<(const exprValue& rhs) const;
 
 
@@ -359,17 +350,12 @@ public:
         //- Return info proxy for printing information to a stream
         InfoProxy<exprValue> info() const { return *this; }
 
-        //- The exprValue as tokens.
-        //  For none : emits pair of brackets.
-        //  For invalid : emits "bad".
-        //
-        //  \param prune suppress the output for none/invalid
-        tokenList tokens(bool prune = false) const;
-
         //- Write the (type-specific) content.
         //  For none : emits pair of brackets.
         //  For invalid : emits "bad".
         //
+        //  Use OTstream for the stream to recover as tokens.
+        //
         //  \param prune suppress the output for none/invalid
         void write(Ostream& os, bool prune = false) const;
 };

src/OpenFOAM/expressions/value/exprValueFieldTag.C

@@ -0,0 +1,307 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023-2024 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 "exprValueFieldTag.H"
+#include "PstreamReduceOps.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::expressions::exprValueFieldTag::empty() const noexcept
+{
+    return
+    (
+        uniformity_ == Foam::Detail::ListPolicy::uniformity::EMPTY
+    );
+}
+
+
+bool Foam::expressions::exprValueFieldTag::is_uniform() const noexcept
+{
+    return
+    (
+        uniformity_ == Foam::Detail::ListPolicy::uniformity::UNIFORM
+    );
+}
+
+
+bool Foam::expressions::exprValueFieldTag::is_nonuniform() const noexcept
+{
+    return
+    (
+        uniformity_ == Foam::Detail::ListPolicy::uniformity::NONUNIFORM
+        // Extra safety for direct reductions?
+        // || uniformity_ == Foam::Detail::ListPolicy::uniformity::MIXED
+    );
+}
+
+
+const Foam::expressions::exprValue&
+Foam::expressions::exprValueFieldTag::value() const noexcept
+{
+    return value_;
+}
+
+
+void Foam::expressions::exprValueFieldTag::set_empty()
+{
+    uniformity_ = Foam::Detail::ListPolicy::uniformity::EMPTY;
+    value_ = Foam::zero{};
+}
+
+
+void Foam::expressions::exprValueFieldTag::set_nouniform()
+{
+    uniformity_ = Foam::Detail::ListPolicy::uniformity::NONUNIFORM;
+    value_ = Foam::zero{};
+}
+
+
+int Foam::expressions::exprValueFieldTag::compare
+(
+    const exprValueFieldTag& rhs
+) const
+{
+    if (uniformity_ != rhs.uniformity_)
+    {
+        // First compare by uniformity
+        return (int(uniformity_) - int(rhs.uniformity_));
+    }
+    if (this == &rhs)
+    {
+        // Identical objects
+        return 0;
+    }
+
+    return value_.compare(rhs.value_);
+}
+
+
+bool Foam::expressions::exprValueFieldTag::equal
+(
+    const exprValueFieldTag& rhs
+) const
+{
+    return (value_ == rhs.value_);
+}
+
+
+void Foam::expressions::exprValueFieldTag::reduce()
+{
+    if (!UPstream::is_parallel())
+    {
+        // Nothing to do
+        return;
+    }
+
+    // Two-stage reduction
+    // ~~~~~~~~~~~~~~~~~~~
+    //
+    // Fields will usually be non-uniform somewhere, so first check with
+    // the cheapest option (bit-wise Allreduce).
+    // Only if they are uniform (with/without empty) do we actually
+    // need to compare values.
+
+    typedef unsigned char bitmask_type;
+
+    bitmask_type shape = static_cast<bitmask_type>(uniformity_);
+
+    // Step 1
+    // ~~~~~~
+    Foam::reduce
+    (
+        shape,
+        bitOrOp<bitmask_type>{},
+        UPstream::msgType(),  // ignored
+        UPstream::worldComm
+    );
+
+    // Step 2
+    // ~~~~~~
+    if
+    (
+        shape == static_cast<bitmask_type>
+        (
+            Foam::Detail::ListPolicy::uniformity::EMPTY
+        )
+    )
+    {
+        // no-op (empty everywhere)
+        value_ = Foam::zero{};
+    }
+    else if
+    (
+        shape == static_cast<bitmask_type>
+        (
+            Foam::Detail::ListPolicy::uniformity::UNIFORM
+        )
+    )
+    {
+        // Ranks are locally uniform (or empty), need to check values too
+        Foam::reduce
+        (
+            *this,
+            exprValueFieldTag::combineOp{},
+            UPstream::msgType(),
+            UPstream::worldComm
+        );
+    }
+    else
+    {
+        // Field is global non-empty and not uniform
+        set_nouniform();
+    }
+}
+
+
+Foam::expressions::exprValueFieldTag
+Foam::expressions::exprValueFieldTag::returnReduce
+(
+    const exprValueFieldTag& tag
+)
+{
+    exprValueFieldTag work(tag);
+    work.reduce();
+    return work;
+}
+
+
+void Foam::expressions::exprValueFieldTag::combine
+(
+    const exprValueFieldTag& b
+)
+{
+    exprValueFieldTag& a = *this;
+
+    if (b.empty())
+    {
+        // no-op
+        return;
+    }
+    else if (a.empty())
+    {
+        a = b;
+    }
+    else if (a.is_uniform() && (!b.is_uniform() || !a.equal(b)))
+    {
+        // Handle two cases:
+        // 1. uniform / non-uniform
+        // 2. uniform / uniform, but with different values
+        a.set_nouniform();
+    }
+
+    // No meaningful value if it is not uniform.
+    // So use zero, but keep the type
+    if (!a.is_uniform())
+    {
+        a.value_ = Foam::zero{};
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+bool Foam::expressions::exprValueFieldTag::
+operator==(const exprValueFieldTag& rhs) const
+{
+    if (uniformity_ != rhs.uniformity_)
+    {
+        // Uniformity must match
+        return false;
+    }
+    else if (this == &rhs)
+    {
+        return true;
+    }
+
+    return (value_ == rhs.value_);
+}
+
+
+bool Foam::expressions::exprValueFieldTag::
+operator<(const exprValueFieldTag& rhs) const
+{
+    return (this->compare(rhs) < 0);
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+void Foam::expressions::exprValueFieldTag::read(Istream& is)
+{
+    label uniformTag;
+
+    is.readBegin("fieldTag");
+
+    is >> uniformTag;
+    uniformity_ = int(uniformTag);
+    value_.read(is);
+
+    is.readEnd("fieldTag");
+}
+
+
+void Foam::expressions::exprValueFieldTag::write(Ostream& os) const
+{
+    os  << token::BEGIN_LIST
+        << label(uniformity_) << token::SPACE;
+    value_.write(os, false);  // No pruning
+    os  << token::END_LIST;
+}
+
+
+void Foam::expressions::exprValueFieldTag::print(Ostream& os) const
+{
+    os  << "{ uniform:"
+        << label(uniformity_)
+        << " type:" << label(value_.typeCode())
+        << " value: " << value_ << " }";
+}
+
+
+Foam::Istream& Foam::operator>>
+(
+    Istream& is,
+    expressions::exprValueFieldTag& tag
+)
+{
+    tag.read(is);
+    return is;
+}
+
+
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const expressions::exprValueFieldTag& tag
+)
+{
+    tag.write(os);
+    return os;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/expressions/value/exprValueFieldTag.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2023 OpenCFD Ltd.
+    Copyright (C) 2023-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,11 +27,11 @@ Class
     Foam::expressions::exprValueFieldTag
 
 Description
-    A polymorphic single-value container for tracking Field content
-    as uniform etc.
+    An expressions::exprValue (polymorphic typed union) with an additional
+    flag for tracking Field content as uniform etc.
 
 SourceFiles
-    exprValueFieldTag.cxx
+    exprValueFieldTag.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -39,7 +39,7 @@ SourceFiles
 #define Foam_expressions_exprValueFieldTag_H
 
 #include "exprValue.H"
-#include "List.H"  // For ListPolicy
+#include "UList.H"  // For ListPolicy
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -54,20 +54,38 @@ namespace expressions
 
 class exprValueFieldTag
 {
-    // Private data
+    // Private Data
 
-        //- Uniformity of field
+        //- Uniformity of field (0: empty, 1: uniform, 2: non-uniform, ...)
+        //  Values as per internal enum Foam::Detail::ListPolicy::uniformity
         int uniformity_{};
 
-        //- Representative (uniform) value for the field
+        //- Representative (uniform) type/value for the field
         expressions::exprValue value_{};
 
+
 public:
 
+    // Generated Methods
+
+        //- Default construct. Uniformity = empty, type = none
+        exprValueFieldTag() = default;
+
+        //- Copy construct
+        exprValueFieldTag(const exprValueFieldTag&) = default;
+
+        //- Copy assignment
+        exprValueFieldTag& operator=(const exprValueFieldTag&) = default;
+
+
     // Constructors
 
-        //- Default construct
-        exprValueFieldTag() = default;
+        //- Construct as uniform with the specified type/value
+        template<class Type>
+        explicit exprValueFieldTag(const Type& val)
+        {
+            set_value(val);
+        }
 
         //- Construct from a range of values
         template<class Type>
@@ -77,6 +95,18 @@ public:
         }
 
 
+    // Factory Methods
+
+        //- Make an empty field tag with the specified type (zero-value)
+        template<class Type>
+        static exprValueFieldTag make_empty()
+        {
+            exprValueFieldTag tag;  // construct empty, no type
+            tag.set_zero<Type>();   // set type and zero value
+            return tag;
+        }
+
+
     // Member Functions
 
         //- True if the uniformity is "empty"
@@ -88,9 +118,18 @@ public:
         //- True if the uniformity is "non-uniform"
         bool is_nonuniform() const noexcept;
 
+        //- Representative (uniform) value for the field
+        const expressions::exprValue& value() const noexcept;
+
+        //- Compare (uniformity, type, value)
+        int compare(const exprValueFieldTag& rhs) const;
+
         //- Test for equality of the values
         bool equal(const exprValueFieldTag& rhs) const;
 
+
+    // Setters
+
         //- Set value and uniformity from range of data
         template<class Type>
         void set(const Type* first, const Type* last)
@@ -104,29 +143,48 @@ public:
             }
             else
             {
+                // Is empty, set zero value
                 value_.set<Type>(Foam::zero{});
             }
         }
 
-        //- Set uniform type and value
+        //- Set as uniform, with specified value
         template<class Type>
-        void set_uniform(const Type& val)
+        void set_value(const Type& val)
         {
             uniformity_ = Foam::Detail::ListPolicy::uniformity::UNIFORM;
             value_.set<Type>(val);
         }
 
-        //- Set as non-uniform
-        void set_nouniform() noexcept;
+        //- Set type and zero value (does not affect uniformity)
+        template<class Type>
+        void set_zero()
+        {
+            value_.set<Type>(Foam::zero{});
+        }
+
+        //- Set as empty with zero value, leave type unchanged
+        void set_empty();
+
+        //- Set as non-uniform with zero value, leave type unchanged
+        void set_nouniform();
+
+
+    // Parallel
+
+        //- Inplace parallel reduction, uses worldComm
+        void reduce();
+
+        //- Perform a reduction on a copy and return the result
+        static exprValueFieldTag returnReduce(const exprValueFieldTag& tag);
 
 
     // Reduction operations
 
-        //- Combine - eg, for global uniformity
+        //- Inplace combine - eg, for global uniformity
         void combine(const exprValueFieldTag& b);
 
-        //- Binary operator to be used by reduce function for detecting
-        //- global uniformity
+        //- Binary combine operator, e.g. for global reduction
         struct combineOp
         {
             exprValueFieldTag operator()
@@ -144,11 +202,32 @@ public:
 
     // IO Operations
 
+        //- Read uniformity label and the value as pair
         void read(Istream& is);
+
+        //- Write uniformity label and the value as pair
         void write(Ostream& os) const;
 
         //- Print description to Ostream
         void print(Ostream& os) const;
+
+
+    // Member Operators
+
+        //- Assign from zero. Changes value but not type
+        void operator=(const Foam::zero) { value_ = Foam::zero{}; }
+
+        //- Compare (uniformity,value) for equality
+        bool operator==(const exprValueFieldTag&) const;
+
+        //- Compare (uniformity,value) for inequality
+        bool operator!=(const exprValueFieldTag& rhs) const
+        {
+            return !(*this == rhs);
+        }
+
+        //- Compare (uniformity,value)
+        bool operator<(const exprValueFieldTag&) const;
 };
 
 } // End namespace expressions

src/OpenFOAM/expressions/value/exprValueFwd.H

@@ -0,0 +1,56 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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.
+
+    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::expressions::exprValue
+
+Description
+    Forward declarations for Foam::expressions::exprValue etc
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_expressions_exprValueFwd_H
+#define Foam_expressions_exprValueFwd_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace expressions
+{
+
+// Forward Declarations
+class exprValue;
+class exprValueFieldTag;
+
+} // End namespace expressions
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/fields/UniformDimensionedFields/UniformDimensionedField.H

@@ -100,8 +100,8 @@ public:
         //- Construct as copy
         UniformDimensionedField(const UniformDimensionedField<Type>&);
 
-        //- Construct from Istream
-        UniformDimensionedField(const IOobject& io);
+        //- Construct from IOobject. Either reads or sets dimensionless zero
+        explicit UniformDimensionedField(const IOobject& io);
 
 
     //- Destructor
@@ -129,7 +129,7 @@ public:
         }
 
         //- Return complete path + object name if the file exists
-        //  either in the case/processor or case otherwise null
+        //- either in the case/processor or case otherwise null
         virtual fileName filePath() const
         {
             return globalFilePath(type());

src/OpenFOAM/matrices/lduMatrix/solvers/GMRES/gmresSolver.C

@@ -1,317 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2010-2015 Hrvoje Jasak, Wikki 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/>.
-
-Description
-    Preconditioned Generalised Minimal Residual solver with
-    run-time selectable preconditioning
-
-\*---------------------------------------------------------------------------*/
-
-#include "gmresSolver.H"
-#include "scalarMatrices.H"
-#include "PrecisionAdaptor.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTypeNameAndDebug(gmresSolver, 0);
-
-    lduMatrix::solver::addsymMatrixConstructorToTable<gmresSolver>
-        addgmresSolverSymMatrixConstructorToTable_;
-
-    lduMatrix::solver::addasymMatrixConstructorToTable<gmresSolver>
-        addgmresSolverAsymMatrixConstructorToTable_;
-
-}
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-void Foam::gmresSolver::givensRotation
-(
-    const solveScalar& h,
-    const solveScalar& beta,
-    solveScalar& c,
-    solveScalar& s
-) const
-{
-    if (beta == 0)
-    {
-        c = 1;
-        s = 0;
-    }
-    else if (mag(beta) > mag(h))
-    {
-        scalar tau = -h/beta;
-        s = 1.0/Foam::sqrt(1.0 + sqr(tau));
-        c = s*tau;
-    }
-    else
-    {
-        scalar tau = -beta/h;
-        c = 1.0/Foam::sqrt(1.0 + sqr(tau));
-        s = c*tau;
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-//- Construct from matrix and solver data stream
-Foam::gmresSolver::gmresSolver
-(
-    const word& fieldName,
-    const lduMatrix& matrix,
-    const FieldField<Field, scalar>& coupleBouCoeffs,
-    const FieldField<Field, scalar>& coupleIntCoeffs,
-    const lduInterfaceFieldPtrsList& interfaces,
-    const dictionary& dict
-)
-:
-    lduMatrix::solver
-    (
-        fieldName,
-        matrix,
-        coupleBouCoeffs,
-        coupleIntCoeffs,
-        interfaces,
-        dict
-    ),
-    preconPtr_
-    (
-        lduMatrix::preconditioner::New
-        (
-            *this,
-            dict
-        )
-    ),
-    nDirs_(dict.getLabel("nDirections"))
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-Foam::solverPerformance Foam::gmresSolver::scalarSolve
-(
-    solveScalarField& x,
-    const solveScalarField& b,
-    const direction cmpt
-) const
-{
-    // Initialise the solverPerformance object to track solver performance
-    solverPerformance solverPerf
-    (
-        preconPtr_->type() + typeName, fieldName()
-    );
-
-    solveScalarField wA(x.size());
-    solveScalarField rA(x.size());
-
-    // Calculate initial residual
-    matrix_.Amul(wA, x, interfaceBouCoeffs_, interfaces_, cmpt);
-
-    // Use rA as scratch space when calculating the normalisation factor
-    solveScalar normFactor = this->normFactor(x, b, wA, rA);
-
-    if (lduMatrix::debug >= 2)
-    {
-        Info<< "   Normalisation factor = " << normFactor << endl;
-    }
-
-    // Calculate residual
-    forAll (rA, i)
-    {
-        rA[i] = b[i] - wA[i];
-    }
-
-    solverPerf.initialResidual() =
-        gSumMag(rA, matrix().mesh().comm())
-       /normFactor;
-    solverPerf.finalResidual() = solverPerf.initialResidual();
-
-    // Note: GMRES cannot be forced to do minIter sweeps
-    // if the residual is zero, due to algorithmic reasons
-    // HJ, 22/Aug/2012
-    if (!solverPerf.checkConvergence(tolerance_, relTol_, log_))
-    {
-        typedef SquareMatrix<solveScalar> solveScalarSquareMatrix;
-
-        // Create the Hesenberg matrix
-        solveScalarSquareMatrix H(nDirs_, 0);
-
-        // Create y and b for Hessenberg matrix
-        solveScalarField yh(nDirs_, 0);
-        solveScalarField bh(nDirs_ + 1, 0);
-
-        // Givens rotation vectors
-        solveScalarField c(nDirs_, 0);
-        solveScalarField s(nDirs_, 0);
-
-        // Allocate Krylov space vectors
-        FieldField<Field, solveScalar> V(nDirs_ + 1);
-
-        forAll (V, i)
-        {
-            V.set(i, new solveScalarField(x.size(), 0));
-        }
-
-        do
-        {
-            // Execute preconditioning
-            preconPtr_->precondition(wA, rA, cmpt);
-
-            // Calculate beta and scale first vector
-            solveScalar beta = Foam::sqrt(gSumSqr(wA, matrix().mesh().comm()));
-
-            // Set initial rhs and bh[0] = beta
-            bh = 0;
-            bh[0] = beta;
-
-            for (label i = 0; i < nDirs_; ++i)
-            {
-                // Set search direction
-                V[i] = wA;
-                V[i] /= beta;
-
-                // Arnoldi's method
-                matrix_.Amul(rA, V[i], interfaceBouCoeffs_, interfaces_, cmpt);
-
-                // Execute preconditioning
-                preconPtr_->precondition(wA, rA, cmpt);
-
-                for (label j = 0; j <= i; ++j)
-                {
-                    beta = gSumProd(wA, V[j], matrix().mesh().comm());
-
-                    H[j][i] = beta;
-
-                    forAll (wA, wI)
-                    {
-                        wA[wI] -= beta*V[j][wI];
-                    }
-                }
-
-                beta = Foam::sqrt(gSumSqr(wA, matrix().mesh().comm()));
-
-                // Apply previous Givens rotations to new column of H.
-                for (label j = 0; j < i; ++j)
-                {
-                    const solveScalar Hji = H[j][i];
-                    H[j][i] = c[j]*Hji - s[j]*H[j + 1][i];
-                    H[j + 1][i] = s[j]*Hji + c[j]*H[j + 1][i];
-                }
-
-                // Apply Givens rotation to current row.
-                givensRotation(H[i][i], beta, c[i], s[i]);
-
-                const solveScalar bhi = bh[i];
-                bh[i] = c[i]*bhi - s[i]*bh[i + 1];
-                bh[i + 1] = s[i]*bhi + c[i]*bh[i + 1];
-                H[i][i] = c[i]*H[i][i] - s[i]*beta;
-            }
-
-            // Back substitute to solve Hy = b
-            for (label i = nDirs_ - 1; i >= 0; i--)
-            {
-                solveScalar sum = bh[i];
-
-                for (label j = i + 1; j < nDirs_; ++j)
-                {
-                    sum -= H[i][j]*yh[j];
-                }
-
-                yh[i] = sum/H[i][i];
-            }
-
-            // Update solution
-            for (label i = 0; i < nDirs_; ++i)
-            {
-                const solveScalarField& Vi = V[i];
-                const solveScalar& yi = yh[i];
-
-                forAll (x, psiI)
-                {
-                    x[psiI] += yi*Vi[psiI];
-                }
-            }
-
-            // Re-calculate the residual
-            matrix_.Amul(wA, x, interfaceBouCoeffs_, interfaces_, cmpt);
-
-            forAll (rA, raI)
-            {
-                rA[raI] = b[raI] - wA[raI];
-            }
-
-            solverPerf.finalResidual() =
-                gSumMag(rA, matrix().mesh().comm())
-               /normFactor;
-            solverPerf.nIterations()++;
-        } while
-        (
-            (
-              ++solverPerf.nIterations() < maxIter_
-            && !solverPerf.checkConvergence(tolerance_, relTol_, log_)
-            )
-         || solverPerf.nIterations() < minIter_
-        );
-    }
-
-
-    if (preconPtr_)
-    {
-        preconPtr_->setFinished(solverPerf);
-    }
-
-    matrix().setResidualField
-    (
-        ConstPrecisionAdaptor<scalar, solveScalar>(rA)(),
-        fieldName_,
-        false
-    );
-
-    return solverPerf;
-}
-
-
-Foam::solverPerformance Foam::gmresSolver::solve
-(
-    scalarField& psi_s,
-    const scalarField& source,
-    const direction cmpt
-) const
-{
-    PrecisionAdaptor<solveScalar, scalar> tpsi(psi_s);
-    return scalarSolve
-    (
-        tpsi.ref(),
-        ConstPrecisionAdaptor<solveScalar, scalar>(source)(),
-        cmpt
-    );
-}
-
-
-// ************************************************************************* //

src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2023 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -95,9 +95,7 @@ Foam::polyPatch::polyPatch
         bm.mesh().points()
     ),
     start_(start),
-    boundaryMesh_(bm),
-    faceCellsPtr_(nullptr),
-    mePtr_(nullptr)
+    boundaryMesh_(bm)
 {
     if (constraintType(patchType))
     {
@@ -124,9 +122,7 @@ Foam::polyPatch::polyPatch
         bm.mesh().points()
     ),
     start_(start),
-    boundaryMesh_(bm),
-    faceCellsPtr_(nullptr),
-    mePtr_(nullptr)
+    boundaryMesh_(bm)
 {}
 
 
@@ -151,9 +147,7 @@ Foam::polyPatch::polyPatch
         bm.mesh().points()
     ),
     start_(dict.get<label>("startFace")),
-    boundaryMesh_(bm),
-    faceCellsPtr_(nullptr),
-    mePtr_(nullptr)
+    boundaryMesh_(bm)
 {
     if (constraintType(patchType))
     {
@@ -180,9 +174,7 @@ Foam::polyPatch::polyPatch
         bm.mesh().points()
     ),
     start_(pp.start()),
-    boundaryMesh_(bm),
-    faceCellsPtr_(nullptr),
-    mePtr_(nullptr)
+    boundaryMesh_(bm)
 {}
 
 
@@ -207,9 +199,7 @@ Foam::polyPatch::polyPatch
         bm.mesh().points()
     ),
     start_(newStart),
-    boundaryMesh_(bm),
-    faceCellsPtr_(nullptr),
-    mePtr_(nullptr)
+    boundaryMesh_(bm)
 {}
 
 
@@ -234,9 +224,7 @@ Foam::polyPatch::polyPatch
         bm.mesh().points()
     ),
     start_(newStart),
-    boundaryMesh_(bm),
-    faceCellsPtr_(nullptr),
-    mePtr_(nullptr)
+    boundaryMesh_(bm)
 {}
 
 
@@ -245,9 +233,7 @@ Foam::polyPatch::polyPatch(const polyPatch& p)
     patchIdentifier(p),
     primitivePatch(p),
     start_(p.start_),
-    boundaryMesh_(p.boundaryMesh_),
-    faceCellsPtr_(nullptr),
-    mePtr_(nullptr)
+    boundaryMesh_(p.boundaryMesh_)
 {}
 
 
@@ -350,25 +336,56 @@ Foam::tmp<Foam::vectorField> Foam::polyPatch::faceCellCentres() const
 }
 
 
-Foam::tmp<Foam::scalarField> Foam::polyPatch::areaFraction() const
+Foam::tmp<Foam::scalarField> Foam::polyPatch::areaFraction
+(
+    const pointField& points
+) const
 {
     auto tfraction = tmp<scalarField>::New(size());
     auto& fraction = tfraction.ref();
 
     const vectorField::subField faceAreas = this->faceAreas();
-    const pointField& points = this->points();
 
     forAll(*this, facei)
     {
-        const face& curFace = this->operator[](facei);
-        fraction[facei] =
-            mag(faceAreas[facei])/(curFace.mag(points) + ROOTVSMALL);
+        const face& f = this->operator[](facei);
+        fraction[facei] = faceAreas[facei].mag()/(f.mag(points) + ROOTVSMALL);
     }
 
     return tfraction;
 }
 
 
+Foam::tmp<Foam::scalarField> Foam::polyPatch::areaFraction() const
+{
+    if (areaFractionPtr_)
+    {
+        return tmp<scalarField>(*areaFractionPtr_);
+    }
+    return nullptr;
+}
+
+
+void Foam::polyPatch::areaFraction(const scalar fraction)
+{
+    areaFractionPtr_ = std::make_unique<scalarField>(size(), fraction);
+}
+
+
+void Foam::polyPatch::areaFraction(const tmp<scalarField>& fraction)
+{
+    if (fraction)
+    {
+        // Steal or clone
+        areaFractionPtr_.reset(fraction.ptr());
+    }
+    else
+    {
+        areaFractionPtr_.reset(nullptr);
+    }
+}
+
+
 const Foam::labelUList& Foam::polyPatch::faceCells() const
 {
     if (!faceCellsPtr_)
@@ -413,6 +430,7 @@ void Foam::polyPatch::clearAddressing()
     primitivePatch::clearPatchMeshAddr();
     faceCellsPtr_.reset(nullptr);
     mePtr_.reset(nullptr);
+    areaFractionPtr_.reset(nullptr);
 }
 
 

src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2015-2023 OpenCFD Ltd.
+    Copyright (C) 2015-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -88,6 +88,9 @@ class polyPatch
         //- Demand-driven: global edge addressing
         mutable std::unique_ptr<labelList> mePtr_;
 
+        //- Cached area fraction
+        std::unique_ptr<scalarField> areaFractionPtr_;
+
 
 protected:
 
@@ -441,10 +444,20 @@ public:
             //- Return face cell centres
             tmp<vectorField> faceCellCentres() const;
 
-            //- Return the area fraction as the ratio of the stored face area
-            //- and the area given by the face points
+            //- Calculate the area fraction as the ratio of the stored face
+            //- area and the area given by the face points.
+            tmp<scalarField> areaFraction(const pointField& points) const;
+
+            //- Return the cached area fraction.
+            //- Usually only set for the non-overlap patches on ACMI.
             tmp<scalarField> areaFraction() const;
 
+            //- Set uniform cached area fraction
+            void areaFraction(const scalar fraction);
+
+            //- Set cached area fraction (non-uniform)
+            void areaFraction(const tmp<scalarField>& fraction);
+
 
         // Addressing into mesh
 

src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2020,2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,7 +28,7 @@ Class
     Foam::tetDecomposer
 
 Description
-    Decomposes polyMesh into tets or pyramids.
+    Decomposes polyMesh into tets (or pyramids)
 
     Cells neighbouring decomposed cells are not decomposed themselves
     so will be polyhedral.
@@ -46,6 +46,7 @@ SourceFiles
 #include "boolList.H"
 #include "typeInfo.H"
 #include "Enum.H"
+#include "faceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -72,7 +73,8 @@ public:
             FACE_CENTRE_TRIS, //- Faces decomposed into triangles
                               //  using face-centre
             FACE_DIAG_TRIS,   //- Faces decomposed into triangles diagonally
-            PYRAMID           //- Faces not decomposed (generates pyramids)
+            PYRAMID,          //- Faces not decomposed (generates pyramids)
+            FACE_DIAG_QUADS   //- Like FACE_DIAG_TRIS but does quads mainly
         };
         static const Enum<decompositionType> decompositionTypeNames;
 
@@ -89,7 +91,6 @@ private:
         //- From face to tet point
         labelList faceToPoint_;
 
-
         // Per face, per point (faceCentre) or triangle (faceDiag)
         // the added tet on the owner side
         labelListList faceOwnerCells_;
@@ -119,6 +120,7 @@ private:
         (
             polyTopoChange& meshMod,
             const face& f,
+            const label facei,
             const label own,
             const label nei,
             const label masterPointID,
@@ -132,6 +134,46 @@ private:
         //- Work out triangle index given the starting vertex in the face
         label triIndex(const label facei, const label fp) const;
 
+        //- Calculate triangulation of boundary faces
+        void splitBoundaryFaces
+        (
+            List<faceList>& boundaryQuads,
+            List<faceList>& boundaryTris
+        ) const;
+
+        //- Correct coupled faces to match up
+        void relativeIndicesToFace
+        (
+            const bool doFlip,
+            const face& meshFace,
+            const faceList& indexLists,
+            faceList& faces
+        ) const;
+
+        //- Calculate triangulation of any face
+        void splitFace
+        (
+            const List<faceList>& boundaryQuads,
+            const List<faceList>& boundaryTris,
+            const label facei,
+            const label patchi,
+            label& quadi,
+            faceList& quadFaces,
+            label& trii,
+            faceList& triFaces
+        ) const;
+
+        void splitFace
+        (
+            const List<faceList>& boundaryQuads,
+            const List<faceList>& boundaryTris,
+            const label facei,
+            const label patchi,
+            faceList& quadFaces,
+            faceList& triFaces,
+            faceList& subFaces
+        ) const;
+
         //- No copy construct
         tetDecomposer(const tetDecomposer&) = delete;
 
@@ -188,7 +230,7 @@ public:
         // Edit
 
             //- Insert all changes into meshMod to convert the polyMesh into
-            //  tets.
+            //  subshapes (tets/prisms)
             void setRefinement
             (
                 const decompositionType decomposeType,
@@ -196,6 +238,17 @@ public:
                 polyTopoChange& meshMod
             );
 
+            //- Insert all changes into meshMod to convert the polyMesh into
+            //- subshapes (tets/prisms). Explicit control over which faces
+            //- get decomposed. Can be used e.g. to not split triangles.
+            void setRefinement
+            (
+                const decompositionType decomposeType,
+                const bitSet& decomposeCell,
+                const bitSet& decomposeFace,
+                polyTopoChange& meshMod
+            );
+
             //- Force recalculation of locally stored data on topological change
             void updateMesh(const mapPolyMesh&);
 };

src/finiteArea/fields/faPatchFields/constraint/empty/emptyFaPatchField.C

@@ -137,4 +137,14 @@ void Foam::emptyFaPatchField<Type>::updateCoeffs()
 }
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::emptyFaPatchField<Type>::write(Ostream& os) const
+{
+    faPatchField<Type>::write(os);
+    // Never write "value"
+}
+
+
 // ************************************************************************* //

src/finiteArea/fields/faPatchFields/constraint/empty/emptyFaPatchField.H

@@ -186,6 +186,12 @@ public:
             {
                 return tmp<Field<Type>>::New();
             }
+
+
+    // Member Functions
+
+        //- Write without "value" entry!
+        virtual void write(Ostream&) const;
 };
 
 

src/finiteArea/fields/faePatchFields/basic/sliced/slicedFaePatchField.C

@@ -131,4 +131,14 @@ Foam::slicedFaePatchField<Type>::~slicedFaePatchField()
 }
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::slicedFaePatchField<Type>::write(Ostream& os) const
+{
+    faePatchField<Type>::write(os);
+    faePatchField<Type>::writeValueEntry(os);
+}
+
+
 // ************************************************************************* //

src/finiteArea/fields/faePatchFields/basic/sliced/slicedFaePatchField.H

@@ -141,6 +141,9 @@ public:
         //- True: this patch field fixes a value.
         virtual bool fixesValue() const { return true; }
 
+        //- Write includes "value" entry
+        virtual void write(Ostream&) const;
+
 
     // Member Operators
 

src/finiteArea/fields/faePatchFields/constraint/empty/emptyFaePatchField.C

@@ -108,4 +108,14 @@ Foam::emptyFaePatchField<Type>::emptyFaePatchField
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::emptyFaePatchField<Type>::write(Ostream& os) const
+{
+    faePatchField<Type>::write(os);
+    // Never write "value"
+}
+
+
 // ************************************************************************* //

src/finiteArea/fields/faePatchFields/constraint/empty/emptyFaePatchField.H

@@ -124,9 +124,9 @@ public:
     virtual ~emptyFaePatchField() = default;
 
 
-    // Member functions
+    // Member Functions
 
-        // Mapping functions
+        // Mapping Functions
 
             //- Map (and resize as needed) from self given a mapping object
             virtual void autoMap
@@ -142,6 +142,12 @@ public:
                 const labelList&
             )
             {}
+
+
+    // Member Functions
+
+        //- Write without "value" entry!
+        virtual void write(Ostream&) const;
 };
 
 

src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaGrad.C

@@ -105,12 +105,13 @@ leastSquaresFaGrad<Type>::calcGrad
     {
         const faPatchField<Type>& bf = vsf.boundaryField()[patchi];
 
-        const Field<Type>& vsfp =
-        (
-            bf.coupled()
-          ? bf.patchNeighbourField().cref()
-          : const_cast<faPatchField<Type>&>(bf)
-        );
+        tmp<Field<Type>> tvsfp(bf);
+
+        if (bf.coupled())
+        {
+            tvsfp = bf.patchNeighbourField();
+        }
+        const auto& vsfp = tvsfp();
 
         const faePatchVectorField& ownLsp = ownLs.boundaryField()[patchi];
         const labelUList& edgeFaces =

src/finiteArea/interpolation/edgeInterpolation/edgeInterpolationScheme/edgeInterpolationScheme.C

@@ -209,8 +209,8 @@ Foam::edgeInterpolationScheme<Type>::interpolate
             label size = vf.boundaryField()[pi].patch().size();
             label start = vf.boundaryField()[pi].patch().start();
 
-            Field<Type> pOwnVf = vf.boundaryField()[pi].patchInternalField();
-            Field<Type> pNgbVf = vf.boundaryField()[pi].patchNeighbourField();
+            Field<Type> pOwnVf(vf.boundaryField()[pi].patchInternalField());
+            Field<Type> pNgbVf(vf.boundaryField()[pi].patchNeighbourField());
 
             Field<Type>& pSf = sf.boundaryFieldRef()[pi];
 

src/finiteVolume/cfdTools/general/MRF/MRFZoneList.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2017 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2021-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,6 +28,7 @@ License
 
 #include "MRFZoneList.H"
 #include "volFields.H"
+#include "surfaceFields.H"
 #include "fixedValueFvsPatchFields.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -216,6 +217,26 @@ Foam::tmp<Foam::volVectorField> Foam::MRFZoneList::DDt
 }
 
 
+Foam::tmp<Foam::surfaceScalarField> Foam::MRFZoneList::phi() const
+{
+    auto tphi = surfaceScalarField::New
+    (
+        "phiMRF",
+        IOobject::NO_REGISTER,
+        mesh_,
+        dimensionedScalar(dimVelocity*dimArea, Zero)
+    );
+    auto& phi = tphi.ref();
+
+    for (const auto& mrf : *this)
+    {
+        mrf.makeAbsolute(phi);
+    }
+
+    return tphi;
+}
+
+
 void Foam::MRFZoneList::makeRelative(volVectorField& U) const
 {
     for (const auto& mrf: *this)

src/finiteVolume/cfdTools/general/MRF/MRFZoneList.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2018 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2021-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -127,6 +127,9 @@ public:
             const volVectorField& U
         ) const;
 
+        //- Return the MRF absolute flux
+        tmp<surfaceScalarField> phi() const;
+
         //- Make the given absolute velocity relative within the MRF region
         void makeRelative(volVectorField& U) const;
 

src/finiteVolume/expressions/PatchFunction1/lookup/Lookup.C

@@ -53,7 +53,7 @@ Foam::Function1Types::Lookup<Type>::Lookup(const Lookup<Type>& rhs)
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
-inline Type Foam::Function1Types::Lookup<Type>::value(const scalar t) const
+Type Foam::Function1Types::Lookup<Type>::value(const scalar t) const
 {
     const objectRegistry& db = function1Base::obr();
     const auto& obj =
@@ -64,7 +64,7 @@ inline Type Foam::Function1Types::Lookup<Type>::value(const scalar t) const
 
 
 template<class Type>
-inline Type Foam::Function1Types::Lookup<Type>::integrate
+Type Foam::Function1Types::Lookup<Type>::integrate
 (
     const scalar t1,
     const scalar t2

src/finiteVolume/fields/fvPatchFields/constraint/empty/emptyFvPatchField.C

@@ -123,4 +123,12 @@ void Foam::emptyFvPatchField<Type>::updateCoeffs()
 }
 
 
+template<class Type>
+void Foam::emptyFvPatchField<Type>::write(Ostream& os) const
+{
+    fvPatchField<Type>::write(os);
+    // Never write "value"
+}
+
+
 // ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/constraint/empty/emptyFvPatchField.H

@@ -42,7 +42,7 @@ Usage
     \verbatim
     <patchName>
     {
-        type            empty;
+        type    empty;
     }
     \endverbatim
 
@@ -194,6 +194,12 @@ public:
             {
                 return tmp<Field<Type>>::New();
             }
+
+
+    // Member Functions
+
+        //- Write without "value" entry!
+        virtual void write(Ostream&) const;
 };
 
 

src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.C

@@ -317,16 +317,10 @@ void Foam::activePressureForceBaffleVelocityFvPatchVectorField::updateCoeffs()
             Info<< "Open fraction = " << openFraction_ << endl;
         }
 
-        scalar areaFraction = 0.0;
-
-        if (opening_)
-        {
-            areaFraction = openFraction_;
-        }
-        else
-        {
-            areaFraction = 1 - openFraction_;
-        }
+        const scalar areaFraction =
+        (
+            opening_ ? openFraction_ : (1 - openFraction_)
+        );
 
         if (patch().boundaryMesh().mesh().moving())
         {
@@ -342,6 +336,7 @@ void Foam::activePressureForceBaffleVelocityFvPatchVectorField::updateCoeffs()
             ).neighbFvPatch().Sf();
         }
 
+
         // Update this wall patch
         vectorField::subField Sfw = patch().patch().faceAreas();
         vectorField newSfw((1 - areaFraction)*initWallSf_);
@@ -350,16 +345,20 @@ void Foam::activePressureForceBaffleVelocityFvPatchVectorField::updateCoeffs()
             Sfw[facei] = newSfw[facei];
         }
         const_cast<scalarField&>(patch().magSf()) = mag(patch().Sf());
+        // Cache fraction
+        const_cast<polyPatch&>(patch().patch()).areaFraction(1-areaFraction);
+
 
         // Update owner side of cyclic
         const_cast<vectorField&>(cyclicPatch.Sf()) = areaFraction*initCyclicSf_;
-
         const_cast<scalarField&>(cyclicPatch.magSf()) = mag(cyclicPatch.Sf());
+        const_cast<polyPatch&>(cyclicPatch.patch()).areaFraction(areaFraction);
+
 
         // Update neighbour side of cyclic
         const_cast<vectorField&>(nbrPatch.Sf()) = areaFraction*nbrCyclicSf_;
-
         const_cast<scalarField&>(nbrPatch.magSf()) = mag(nbrPatch.Sf());
+        const_cast<polyPatch&>(nbrPatch.patch()).areaFraction(areaFraction);
 
         curTimeIndex_ = this->db().time().timeIndex();
     }

src/finiteVolume/fields/fvsPatchFields/basic/calculated/calculatedFvsPatchField.H

@@ -43,8 +43,8 @@ Usage
     \verbatim
     <patchName>
     {
-        type            calculated;
-        value           uniform (0 0 0);    // Required value entry
+        type    calculated;
+        value   uniform (0 0 0);    // Required value entry
     }
     \endverbatim
 
@@ -88,6 +88,7 @@ public:
         );
 
         //- Construct from patch, internal field and dictionary
+        //- The "value" entry is MUST_READ.
         calculatedFvsPatchField
         (
             const fvPatch&,

src/finiteVolume/fields/fvsPatchFields/basic/coupled/coupledFvsPatchField.H

@@ -123,7 +123,7 @@ public:
 
     // Member Functions
 
-        //- True if this patch field is derived from coupledFvsPatchField
+        //- True: this patch field is derived from coupledFvsPatchField
         virtual bool coupled() const
         {
             return true;

src/finiteVolume/fields/fvsPatchFields/basic/fixedValue/fixedValueFvsPatchField.H

@@ -43,8 +43,8 @@ Usage
     \verbatim
     <patchName>
     {
-        type            fixedValue;
-        value           uniform 0;  // Example for scalar field usage
+        type    fixedValue;
+        value   uniform 0;  // Example for scalar field usage
     }
     \endverbatim
 
@@ -151,7 +151,7 @@ public:
         virtual void write(Ostream&) const;
 
 
-        // Evaluation functions
+        // Evaluation Functions
 
             //- Return the matrix diagonal coefficients corresponding to the
             //  evaluation of the value of this patchField with given weights

src/finiteVolume/fields/fvsPatchFields/constraint/cyclic/cyclicFvsPatchField.H

@@ -88,7 +88,7 @@ public:
         );
 
         //- Construct from patch, internal field and dictionary.
-        //- The "value" entry is NO_READ.
+        //- The "value" entry is MUST_READ.
         cyclicFvsPatchField
         (
             const fvPatch&,

src/finiteVolume/fields/fvsPatchFields/constraint/cyclicSlip/cyclicSlipFvsPatchField.H

@@ -54,7 +54,6 @@ class cyclicSlipFvsPatchField
 :
     public cyclicFvsPatchField<Type>
 {
-
 public:
 
     //- Runtime type information

src/finiteVolume/fields/fvsPatchFields/constraint/empty/emptyFvsPatchField.C

@@ -110,4 +110,14 @@ Foam::emptyFvsPatchField<Type>::emptyFvsPatchField
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::emptyFvsPatchField<Type>::write(Ostream& os) const
+{
+    fvsPatchField<Type>::write(os);
+    // Never write "value"
+}
+
+
 // ************************************************************************* //

src/finiteVolume/fields/fvsPatchFields/constraint/empty/emptyFvsPatchField.H

@@ -37,7 +37,7 @@ Usage
     \verbatim
     <patchName>
     {
-        type            empty;
+        type    empty;
     }
     \endverbatim
 
@@ -146,6 +146,12 @@ public:
                 const labelList&
             )
             {}
+
+
+    // Member Functions
+
+        //- Write without "value" entry!
+        virtual void write(Ostream&) const;
 };
 
 

src/finiteVolume/fields/fvsPatchFields/constraint/nonuniformTransformCyclic/nonuniformTransformCyclicFvsPatchField.H

@@ -54,7 +54,6 @@ class nonuniformTransformCyclicFvsPatchField
 :
     public cyclicFvsPatchField<Type>
 {
-
 public:
 
     //- Runtime type information

src/finiteVolume/fields/fvsPatchFields/constraint/symmetry/symmetryFvsPatchField.C

@@ -49,7 +49,7 @@ Foam::symmetryFvsPatchField<Type>::symmetryFvsPatchField
     const dictionary& dict
 )
 :
-    fvsPatchField<Type>(p, iF, dict, IOobjectOption::NO_READ)
+    fvsPatchField<Type>(p, iF, dict, IOobjectOption::MUST_READ)
 {
     if (!isType<symmetryFvPatch>(p))
     {
@@ -105,4 +105,14 @@ Foam::symmetryFvsPatchField<Type>::symmetryFvsPatchField
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::symmetryFvsPatchField<Type>::write(Ostream& os) const
+{
+    fvsPatchField<Type>::write(os);
+    fvsPatchField<Type>::writeValueEntry(os);
+}
+
+
 // ************************************************************************* //

src/finiteVolume/fields/fvsPatchFields/constraint/symmetry/symmetryFvsPatchField.H

@@ -30,14 +30,15 @@ Class
 Description
     This boundary condition enforces a symmetry constraint
 
-    The "value" entry is NO_READ, NO_WRITE.
+    The "value" entry is MUST_READ.
 
 Usage
     Example of the boundary condition specification:
     \verbatim
     <patchName>
     {
-        type            symmetry;
+        type    symmetry;
+        value   ...;
     }
     \endverbatim
 
@@ -82,7 +83,7 @@ public:
         );
 
         //- Construct from patch, internal field and dictionary
-        //- The "value" entry is NO_READ.
+        //- The "value" entry is MUST_READ.
         symmetryFvsPatchField
         (
             const fvPatch&,
@@ -126,6 +127,12 @@ public:
         {
             return fvsPatchField<Type>::Clone(*this, iF);
         }
+
+
+    // Member Functions
+
+        //- Write includes "value" entry
+        virtual void write(Ostream&) const;
 };
 
 

src/finiteVolume/fields/fvsPatchFields/constraint/symmetryPlane/symmetryPlaneFvsPatchField.C

@@ -49,7 +49,7 @@ Foam::symmetryPlaneFvsPatchField<Type>::symmetryPlaneFvsPatchField
     const dictionary& dict
 )
 :
-    fvsPatchField<Type>(p, iF, dict, IOobjectOption::NO_READ)
+    fvsPatchField<Type>(p, iF, dict, IOobjectOption::MUST_READ)
 {
     if (!isType<symmetryPlaneFvPatch>(p))
     {
@@ -105,4 +105,14 @@ Foam::symmetryPlaneFvsPatchField<Type>::symmetryPlaneFvsPatchField
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::symmetryPlaneFvsPatchField<Type>::write(Ostream& os) const
+{
+    fvsPatchField<Type>::write(os);
+    fvsPatchField<Type>::writeValueEntry(os);
+}
+
+
 // ************************************************************************* //

src/finiteVolume/fields/fvsPatchFields/constraint/symmetryPlane/symmetryPlaneFvsPatchField.H

@@ -29,14 +29,15 @@ Class
 Description
     This boundary condition enforces a symmetryPlane constraint
 
-    The "value" entry is NO_READ, NO_WRITE.
+    The "value" entry is MUST_READ.
 
 Usage
     Example of the boundary condition specification:
     \verbatim
     <patchName>
     {
-        type            symmetryPlane;
+        type    symmetryPlane;
+        value   ...;
     }
     \endverbatim
 
@@ -81,7 +82,7 @@ public:
         );
 
         //- Construct from patch, internal field and dictionary.
-        //- The "value" entry is NO_READ.
+        //- The "value" entry is MUST_READ.
         symmetryPlaneFvsPatchField
         (
             const fvPatch&,
@@ -125,6 +126,12 @@ public:
         {
             return fvsPatchField<Type>::Clone(*this, iF);
         }
+
+
+    // Member Functions
+
+        //- Write includes "value" entry
+        virtual void write(Ostream&) const;
 };
 
 

src/finiteVolume/fields/fvsPatchFields/constraint/wedge/wedgeFvsPatchField.C

@@ -49,7 +49,7 @@ Foam::wedgeFvsPatchField<Type>::wedgeFvsPatchField
     const dictionary& dict
 )
 :
-    fvsPatchField<Type>(p, iF, dict, IOobjectOption::NO_READ)
+    fvsPatchField<Type>(p, iF, dict, IOobjectOption::MUST_READ)
 {
     if (!isType<wedgeFvPatch>(p))
     {
@@ -105,4 +105,14 @@ Foam::wedgeFvsPatchField<Type>::wedgeFvsPatchField
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::wedgeFvsPatchField<Type>::write(Ostream& os) const
+{
+    fvsPatchField<Type>::write(os);
+    fvsPatchField<Type>::writeValueEntry(os);
+}
+
+
 // ************************************************************************* //

src/finiteVolume/fields/fvsPatchFields/constraint/wedge/wedgeFvsPatchField.H

@@ -31,14 +31,15 @@ Description
     This boundary condition is similar to the cyclic condition, except that
     it is applied to 2-D geometries.
 
-    The "value" entry is NO_READ, NO_WRITE.
+    The "value" entry is MUST_READ.
 
 Usage
     Example of the boundary condition specification:
     \verbatim
     <patchName>
     {
-        type            wedge;
+        type    wedge;
+        value   uniform ...;    // Required value entry
     }
     \endverbatim
 
@@ -83,7 +84,7 @@ public:
         );
 
         //- Construct from patch, internal field and dictionary.
-        //- The "value" entry is NO_READ.
+        //- The "value" entry is MUST_READ.
         wedgeFvsPatchField
         (
             const fvPatch&,
@@ -127,6 +128,12 @@ public:
         {
             return fvsPatchField<Type>::Clone(*this, iF);
         }
+
+
+    // Member Functions
+
+        //- Write includes "value" entry
+        virtual void write(Ostream&) const;
 };
 
 

src/finiteVolume/fvMesh/wallDist/wallDistAddressing/wallDistAddressing.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2023 OpenCFD Ltd.
+    Copyright (C) 2023-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -146,29 +146,40 @@ void Foam::wallDistAddressing::correct(volScalarField& y)
     globalWallsPtr_.reset(new globalIndex(nWalls));
     const globalIndex& globalWalls = globalWallsPtr_();
 
-    labelList seedFaces(nWalls);
-    List<wallPointAddressing> seedInfo(nWalls);
+    DynamicList<label> seedFaces(nWalls);
+    DynamicList<wallPointAddressing> seedInfo(nWalls);
 
 
     nWalls = 0;
     for (const label patchi : patchIDs_)
     {
         const auto& fc = C.boundaryField()[patchi];
+        const auto areaFraction(fc.patch().patch().areaFraction());
 
         forAll(fc, patchFacei)
         {
-            seedFaces[nWalls] = fc.patch().start()+patchFacei;
-            seedInfo[nWalls] = wallPointAddressing
+            if
             (
-                fc[patchFacei],
-                gt.encode
+                !areaFraction
+             || (areaFraction()[patchFacei] > 0.5)    // mostly wall
+            )
+            {
+                seedFaces.append(fc.patch().start()+patchFacei);
+                seedInfo.append
                 (
-                    Pstream::myProcNo(),
-                    nWalls,
-                    gt.nullTransformIndex()
-                ),
-                scalar(0.0)
-            );
+                    wallPointAddressing
+                    (
+                        fc[patchFacei],
+                        gt.encode
+                        (
+                            Pstream::myProcNo(),
+                            nWalls,
+                            gt.nullTransformIndex()
+                        ),
+                        scalar(0.0)
+                    )
+                );
+            }
             nWalls++;
         }
     }

src/functionObjects/field/Make/files

@@ -120,6 +120,7 @@ fluxSummary/fluxSummary.C
 mapFields/mapFields.C
 reactionSensitivityAnalysis/reactionsSensitivityAnalysisObjects.C
 DESModelRegions/DESModelRegions.C
+cellDecomposer/cellDecomposer.C
 
 externalCoupled/externalCoupled.C
 externalCoupled/externalCoupledMixed/externalCoupledMixedFvPatchFields.C

src/functionObjects/field/Make/options

@@ -6,6 +6,7 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
@@ -33,6 +34,7 @@ LIB_LIBS = \
     -lmeshTools \
     -ldynamicMesh \
     -lsampling \
+    -lfvOptions \
     -llagrangian \
     -ldistributionModels \
     -lfluidThermophysicalModels \

src/functionObjects/field/cellDecomposer/cellDecomposer.C

@@ -0,0 +1,415 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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.
+
+    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 "cellDecomposer.H"
+#include "addToRunTimeSelectionTable.H"
+#include "polyTopoChange.H"
+#include "mapPolyMesh.H"
+#include "tetDecomposer.H"
+#include "syncTools.H"
+#include "dummyTransform.H"
+#include "ReadFields.H"
+#include "surfaceFields.H"
+#include "PackedBoolList.H"
+#include "fvMeshTools.H"
+#include "cellSetOption.H"
+#include "cellBitSet.H"
+#include "cellSet.H"
+#include "hexMatcher.H"
+#include "prismMatcher.H"
+#include "pyrMatcher.H"
+#include "tetMatcher.H"
+
+#include "OBJstream.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+    defineTypeNameAndDebug(cellDecomposer, 0);
+    addToRunTimeSelectionTable(functionObject, cellDecomposer, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::functionObjects::cellDecomposer::makeMesh
+(
+    const dictionary& dict,
+    const word& regionName
+)
+{
+    Info<< name() << ':' << nl
+        << "    Decomposing cells to region " << regionName << endl;
+
+    Info<< incrIndent;
+
+    const fv::cellSetOption::selectionModeType selectionMode
+    (
+        fv::cellSetOption::selectionModeTypeNames_.get
+        (
+            "selectionMode",
+            dict
+        )
+    );
+
+
+    // Start off from existing mesh
+    polyTopoChange meshMod(mesh_);
+
+    tetDecompPtr_.reset(new Foam::tetDecomposer(mesh_));
+
+    // Default values
+    bitSet decomposeCell(mesh_.nCells());
+    autoPtr<bitSet> decomposeFacePtr;
+    tetDecomposer::decompositionType decompType
+    (
+        tetDecomposer::FACE_CENTRE_TRIS
+    );
+
+
+    switch (selectionMode)
+    {
+        case fv::cellSetOption::smAll:
+        {
+            Info<< indent << "- selecting all cells" << endl;
+
+            decomposeCell = true;
+            break;
+        }
+        case fv::cellSetOption::smGeometric:
+        {
+            Info<< indent << "- selecting cells geometrically" << endl;
+
+            decomposeCell =
+                cellBitSet::select(mesh_, dict.subDict("selection"), true);
+            break;
+        }
+        case fv::cellSetOption::smCellSet:
+        {
+            const word selectionName(dict.get<word>("cellSet"));
+
+            Info<< indent
+                << "- selecting cells using cellSet "
+                << selectionName << endl;
+
+            decomposeCell.set
+            (
+                cellSet(mesh_, selectionName, IOobject::MUST_READ).toc()
+            );
+            break;
+        }
+        case fv::cellSetOption::smCellZone:
+        {
+            wordRes selectionNames;
+            if
+            (
+                !dict.readIfPresent("cellZones", selectionNames)
+             || selectionNames.empty()
+            )
+            {
+                selectionNames.resize(1);
+                dict.readEntry("cellZone", selectionNames.first());
+            }
+
+            Info<< indent
+                << "- selecting cells using cellZones "
+                << flatOutput(selectionNames) << nl;
+
+            const auto& zones = mesh_.cellZones();
+
+            // Also handles groups, multiple zones etc ...
+            labelList zoneIDs = zones.indices(selectionNames);
+
+            if (zoneIDs.empty())
+            {
+                FatalErrorInFunction
+                    << "No matching cellZones: "
+                    << flatOutput(selectionNames) << nl
+                    << "Valid zones : "
+                    << flatOutput(zones.names()) << nl
+                    << "Valid groups: "
+                    << flatOutput(zones.groupNames())
+                    << nl
+                    << exit(FatalError);
+            }
+
+            if (zoneIDs.size() == 1)
+            {
+                decomposeCell.set(zones[zoneIDs.first()]);
+                // TBD: Foam::sort(cells_);
+            }
+            else
+            {
+                decomposeCell.set(zones.selection(zoneIDs).sortedToc());
+            }
+            break;
+        }
+        default:
+        {
+            FatalErrorInFunction
+                << "Unsupported selectionMode "
+                << fv::cellSetOption::selectionModeTypeNames_[selectionMode]
+                << ". Valid selectionMode types are "
+                << fv::cellSetOption::selectionModeTypeNames_
+                << exit(FatalError);
+        }
+    }
+
+    word decompTypeName;
+    if (dict.readIfPresent("decomposeType", decompTypeName))
+    {
+        if (decompTypeName == "polyhedral")
+        {
+            // Automatic selection to generate hex/prism/tet only
+            //  - subset decomposeCell to exclude hex/prism/tet
+            //  - set faces to FACE_DIAG_QUADS
+            // Usually start with cellSetOption::smAll
+            decomposeFacePtr.reset(new bitSet(mesh_.nFaces()));
+            auto& decomposeFace = decomposeFacePtr();
+
+            decompType = tetDecomposer::FACE_DIAG_QUADS;
+            bitSet oldDecomposeCell(decomposeCell);
+            decomposeCell = false;
+
+            // Construct shape recognizers
+            prismMatcher prism;
+
+            for (const label celli : oldDecomposeCell)
+            {
+                if
+                (
+                   !hexMatcher::test(mesh_, celli)
+                && !tetMatcher::test(mesh_, celli)
+                && !pyrMatcher::test(mesh_, celli)
+                && !prism.isA(mesh_, celli)
+                )
+                {
+                    decomposeCell.set(celli);
+                    decomposeFace.set(mesh_.cells()[celli]);
+                }
+            }
+
+            // Sync boundary info
+            syncTools::syncFaceList
+            (
+                mesh_,
+                decomposeFace,
+                orEqOp<unsigned int>()
+            );
+        }
+        else
+        {
+            decompType = tetDecomposer::decompositionTypeNames[decompTypeName];
+        }
+    }
+
+
+    // Insert all changes to create tets
+    if (decomposeFacePtr)
+    {
+        if (debug)
+        {
+            OBJstream os(mesh_.time().path()/"orig_faces.obj");
+            os.write
+            (
+                UIndirectList<face>
+                (
+                    mesh_.faces(),
+                    decomposeFacePtr().sortedToc()
+                )(),
+                mesh_.points(),
+                false
+            );
+            Pout<< "Written " << meshMod.faces().size()
+                << " faces to " << os.name() << endl;
+        }
+
+        tetDecompPtr_().setRefinement
+        (
+            decompType, //Foam::tetDecomposer::FACE_CENTRE_TRIS,
+            decomposeCell,
+            decomposeFacePtr(),
+            meshMod
+        );
+    }
+    else
+    {
+        tetDecompPtr_().setRefinement
+        (
+            decompType, //Foam::tetDecomposer::FACE_CENTRE_TRIS,
+            decomposeCell,
+            meshMod
+        );
+    }
+
+
+    if (debug)
+    {
+        OBJstream os(mesh_.time().path()/"faces.obj");
+        os.write(meshMod.faces(), pointField(meshMod.points()), false);
+        Pout<< "Written " << meshMod.faces().size()
+            << " faces to " << os.name() << endl;
+    }
+
+
+    autoPtr<fvMesh> tetMeshPtr;
+
+    mapPtr_ = meshMod.makeMesh
+    (
+        tetMeshPtr,
+        IOobject
+        (
+            regionName,
+            mesh_.facesInstance(),      // same instance as original mesh
+            mesh_.time(),               //? why same database? TBD.
+            IOobject::READ_IF_PRESENT,  // Read fv* if present
+            IOobject::AUTO_WRITE,
+            IOobject::REGISTER
+        ),
+        mesh_
+    );
+
+
+    //- Update numbering on tet-decomposition engine
+    tetDecompPtr_().updateMesh(mapPtr_());
+
+    Info<< indent << "Writing decomposed mesh to "
+        << tetMeshPtr().objectRegistry::objectRelPath()
+        << endl;
+    tetMeshPtr().write();
+
+    fvMeshTools::createDummyFvMeshFiles(mesh_.time(), regionName, true);
+
+    // Store new mesh on object registry
+    tetMeshPtr.ptr()->polyMesh::store();
+
+    Info<< decrIndent;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::functionObjects::cellDecomposer::cellDecomposer
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    fvMeshFunctionObject(name, runTime, dict)
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::functionObjects::cellDecomposer::read(const dictionary& dict)
+{
+    if (fvMeshFunctionObject::read(dict))
+    {
+        // Generate new tet equivalent mesh. TBD: why meshObject?
+        //meshObjects::tetDecomposition::New(mesh, dict);
+
+        dict_ = dict.optionalSubDict(typeName + "Coeffs");
+        dict_.readEntry("mapRegion", mapRegion_);
+        dict_.readEntry("fields", fieldNames_);
+        makeMesh(dict_, mapRegion_);
+    }
+
+    return true;
+}
+
+
+bool Foam::functionObjects::cellDecomposer::execute()
+{
+    Log << type() << " " << name() << " execute:" << nl;
+
+    bool ok = false;
+
+    if (mesh_.changing())
+    {
+        // Re-do mesh. Currently rebuilds whole mesh. Could move points only
+        // for mesh motion.
+        tetDecompPtr_.clear();
+        mapPtr_.clear();
+        const_cast<Time&>(this->mesh_.time()).erase(mapRegion_);
+        makeMesh(dict_, mapRegion_);
+    }
+
+
+    // Look up
+    ok = mapFieldType<scalar>() || ok;
+    ok = mapFieldType<vector>() || ok;
+    ok = mapFieldType<sphericalTensor>() || ok;
+    ok = mapFieldType<symmTensor>() || ok;
+    ok = mapFieldType<tensor>() || ok;
+
+    if (log)
+    {
+        if (!ok)
+        {
+            Info<< "    none" << nl;
+        }
+
+        Info<< endl;
+    }
+    return true;
+}
+
+
+bool Foam::functionObjects::cellDecomposer::write()
+{
+    Log << type() << " " << name() << " write:" << nl;
+
+    bool ok = false;
+
+    ok = writeFieldType<scalar>() || ok;
+    ok = writeFieldType<vector>() || ok;
+    ok = writeFieldType<sphericalTensor>() || ok;
+    ok = writeFieldType<symmTensor>() || ok;
+    ok = writeFieldType<tensor>() || ok;
+
+    if (log)
+    {
+        if (!ok)
+        {
+            Info<< "    none" << nl;
+        }
+
+        Info<< endl;
+    }
+
+    return true;
+}
+
+
+// ************************************************************************* //

src/functionObjects/field/cellDecomposer/cellDecomposer.H

@@ -0,0 +1,226 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::functionObjects::cellDecomposer
+
+Group
+    grpFieldFunctionObjects
+
+Description
+    Maps input fields from local mesh to a secondary mesh at runtime.
+
+    The secondary mesh gets created on-the-fly by decomposing the current mesh.
+
+
+    Operands:
+    \table
+      Operand      | Type                       | Location
+      input        | {vol,surface}\<Type\>Field <!--
+                                    --> | $FOAM_CASE/\<time\>/\<inpField\>
+      output file  | -                          | -
+      output field | {vol,surface}\<Type\>Field <!--
+                                    --> | $FOAM_CASE/\<time\>/\<outField\>
+    \endtable
+
+    where \c \<Type\>=Scalar/Vector/SphericalTensor/SymmTensor/Tensor.
+
+Usage
+    Minimal example by using \c system/controlDict.functions:
+    \verbatim
+    cellDecomposer1
+    {
+        // Mandatory entries (unmodifiable)
+        type            cellDecomposer;
+        libs            (fieldFunctionObjects);
+
+        // Mandatory (inherited) entries (runtime modifiable)
+        fields          (<field1> <field2> ... <fieldN>);
+        mapRegion       myTetMesh;
+        ...
+        // Mandatory entries
+        // Decompose type: 
+        decomposeType   polyhedral;
+        // Cell set to decompose
+        selectionMode   all;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property     | Description                        | Type | Req'd | Dflt
+      type         | Type name: cellDecomposer          | word |  yes  | -
+      libs         | Library name: fieldFunctionObjects | word |  yes  | -
+      fields       | Names of operand fields            | wordList |  yes  | -
+      mapRegion    | Name of region to map to           | word |  yes  | -
+      decomposeType | How to decompose cells            | word |  yes  | -
+      selectionMode | How to select cells (see fvOption)| word |  yes  | -
+    \endtable
+
+    decomposeType:
+    - faceCentre : decompose cells into tets using face centre and cell centre.
+                   (hex becomes 6*4 tets)
+    - faceDiagonal : decompose cells into tets using face diagonal, similar
+                     to implicit decomposition inside lagrangian tracking.
+                     (hex becomes 6*2 tets)
+    - pyramid : keep faces intact but create (polygonal-base) pyramids using
+                cell centre (hex becomes 6 pyramids)
+    - faceDiagonalQuads : like faceDiagonal but split faces into quads and
+                          triangles instead of just triangles
+    - polyhedral : like faceDiagonalQuads but only decompose non-hex/prism/tet
+                   cells in selected set. Used to convert polyhedral mesh into
+                   'simple' mesh.
+
+    The inherited entries are elaborated in:
+     - \link functionObject.H \endlink
+
+See also
+    - Foam::functionObjects::mapFields
+
+SourceFiles
+    cellDecomposer.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef functionObjects_cellDecomposer_H
+#define functionObjects_cellDecomposer_H
+
+#include "fvMeshFunctionObject.H"
+#include "volFieldsFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class tetDecomposer;
+class mapPolyMesh;
+
+namespace functionObjects
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class cellDecomposer Declaration
+\*---------------------------------------------------------------------------*/
+
+class cellDecomposer
+:
+    public fvMeshFunctionObject
+{
+    // Private Data
+
+        //- Parameter dictionary
+        dictionary dict_;
+
+        //- Name of new mesh
+        word mapRegion_;
+
+        //- List of field names to interpolate
+        wordRes fieldNames_;
+
+        //- Tet decomposer
+        autoPtr<tetDecomposer> tetDecompPtr_;
+
+        //- Map from polyMesh to tet-mesh
+        autoPtr<mapPolyMesh> mapPtr_;
+
+
+    // Private Member Functions
+
+        //- Generate mesh
+        void makeMesh(const dictionary& dict, const word& name);
+
+        //- Helper function to map the \<Type\> fields
+        template<class Type>
+        bool mapFieldType() const;
+
+        //- Helper function to write the \<Type\> fields
+        template<class Type>
+        bool writeFieldType() const;
+
+        template<class Type>
+        tmp<GeometricField<Type, fvPatchField, volMesh>>
+        interpolate
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& vf,
+            const fvMesh& sMesh,
+            const labelUList& patchMap,
+            const labelUList& cellMap,
+            const labelUList& faceMap,
+            const bool allowUnmapped
+        ) const;
+
+
+public:
+
+
+    //- Runtime type information
+    TypeName("cellDecomposer");
+
+
+    // Constructors
+
+        //- Construct from Time and dictionary
+        cellDecomposer
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~cellDecomposer() = default;
+
+
+    // Member Functions
+
+        //- Read the cellDecomposer data
+        virtual bool read(const dictionary& dict);
+
+        //- Execute
+        virtual bool execute();
+
+        //- Write
+        virtual bool write();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "cellDecomposerTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/functionObjects/field/cellDecomposer/cellDecomposerTemplates.C

@@ -0,0 +1,364 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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.
+
+    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 "fvMesh.H"
+#include "emptyFvPatchFields.H"
+#include "directFvPatchFieldMapper.H"
+#include "mapPolyMesh.H"
+//#include "polyPatch.H"
+//#include "lduSchedule.H"
+//#include "meshToMesh.H"
+
+//template<class Type>
+//void Foam::functionObjects::cellDecomposer::evaluateConstraintTypes
+//(
+//    GeometricField<Type, fvPatchField, volMesh>& fld
+//) const
+//{
+//    auto& fldBf = fld.boundaryFieldRef();
+//
+//    const UPstream::commsTypes commsType = UPstream::defaultCommsType;
+//    const label startOfRequests = UPstream::nRequests();
+//
+//    if
+//    (
+//        commsType == UPstream::commsTypes::blocking
+//     || commsType == UPstream::commsTypes::nonBlocking
+//    )
+//    {
+//        forAll(fldBf, patchi)
+//        {
+//            fvPatchField<Type>& tgtField = fldBf[patchi];
+//
+//            if
+//            (
+//                tgtField.type() == tgtField.patch().patch().type()
+//             && polyPatch::constraintType(tgtField.patch().patch().type())
+//            )
+//            {
+//                tgtField.initEvaluate(commsType);
+//            }
+//        }
+//
+//        // Wait for outstanding requests
+//        if (commsType == UPstream::commsTypes::nonBlocking)
+//        {
+//            UPstream::waitRequests(startOfRequests);
+//        }
+//
+//        forAll(fldBf, patchi)
+//        {
+//            fvPatchField<Type>& tgtField = fldBf[patchi];
+//
+//            if
+//            (
+//                tgtField.type() == tgtField.patch().patch().type()
+//             && polyPatch::constraintType(tgtField.patch().patch().type())
+//            )
+//            {
+//                tgtField.evaluate(commsType);
+//            }
+//        }
+//    }
+//    else if (commsType == UPstream::commsTypes::scheduled)
+//    {
+//        const lduSchedule& patchSchedule =
+//            fld.mesh().globalData().patchSchedule();
+//
+//        for (const auto& schedEval : patchSchedule)
+//        {
+//            const label patchi = schedEval.patch;
+//
+//            fvPatchField<Type>& tgtField = fldBf[patchi];
+//
+//            if
+//            (
+//                tgtField.type() == tgtField.patch().patch().type()
+//             && polyPatch::constraintType(tgtField.patch().patch().type())
+//            )
+//            {
+//                if (schedEval.init)
+//                {
+//                    tgtField.initEvaluate(commsType);
+//                }
+//                else
+//                {
+//                    tgtField.evaluate(commsType);
+//                }
+//            }
+//        }
+//    }
+//}
+
+template<class Type>
+Foam::tmp
+<
+    Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>
+>
+Foam::functionObjects::cellDecomposer::interpolate
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vf,
+    const fvMesh& sMesh,
+    const labelUList& patchMap,
+    const labelUList& cellMap,
+    const labelUList& faceMap,
+    const bool allowUnmapped
+) const
+{
+    // 1. Create the complete field with dummy patch fields
+    PtrList<fvPatchField<Type>> patchFields(patchMap.size());
+
+    forAll(patchFields, patchi)
+    {
+        // Set the first one by hand as it corresponds to the
+        // exposed internal faces. Additional interpolation can be put here
+        // as necessary.
+        if (patchMap[patchi] == -1)
+        {
+            patchFields.set
+            (
+                patchi,
+                new emptyFvPatchField<Type>
+                (
+                    sMesh.boundary()[patchi],
+                    fvPatchField<Type>::Internal::null()
+                )
+            );
+        }
+        else
+        {
+            patchFields.set
+            (
+                patchi,
+                fvPatchField<Type>::New
+                (
+                    fvPatchFieldBase::calculatedType(),
+                    sMesh.boundary()[patchi],
+                    fvPatchField<Type>::Internal::null()
+                )
+            );
+        }
+    }
+
+    auto tresult = tmp<GeometricField<Type, fvPatchField, volMesh>>::New
+    (
+        IOobject
+        (
+            "subset"+vf.name(),
+            sMesh.time().timeName(),
+            sMesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        sMesh,
+        vf.dimensions(),
+        Field<Type>(vf.primitiveField(), cellMap),
+        patchFields
+    );
+    auto& result = tresult.ref();
+    result.oriented() = vf.oriented();
+
+
+    // 2. Change the fvPatchFields to the correct type using a mapper
+    //  constructor (with reference to the now correct internal field)
+
+    auto& bf = result.boundaryFieldRef();
+
+    forAll(bf, patchi)
+    {
+        const label basePatchId = patchMap[patchi];
+
+        if (basePatchId != -1)
+        {
+            // Construct addressing
+            const fvPatch& subPatch = sMesh.boundary()[patchi];
+            const fvPatch& basePatch = vf.mesh().boundary()[basePatchId];
+            const label baseStart = basePatch.start();
+            const label baseSize = basePatch.size();
+
+            labelList directAddressing(subPatch.size());
+
+            forAll(directAddressing, i)
+            {
+                const label baseFacei = faceMap[subPatch.start()+i];
+
+                if (baseFacei >= baseStart && baseFacei < baseStart+baseSize)
+                {
+                    directAddressing[i] = baseFacei-baseStart;
+                }
+                else
+                {
+                    // Mapped from internal face. Do what? Leave up to
+                    // fvPatchField
+                    directAddressing[i] = -1;
+                }
+            }
+
+
+            directFvPatchFieldMapper mapper(directAddressing);
+
+            // allowUnmapped : special mode for if we do not want to be
+            // warned for unmapped faces (e.g. from fvMeshDistribute).
+
+            const bool hasUnmapped = mapper.hasUnmapped();
+            if (allowUnmapped)
+            {
+                mapper.hasUnmapped() = false;
+            }
+
+            bf.set
+            (
+                patchi,
+                fvPatchField<Type>::New
+                (
+                    vf.boundaryField()[basePatchId],
+                    subPatch,
+                    result.internalField(),
+                    mapper
+                )
+            );
+
+            if (allowUnmapped && hasUnmapped)
+            {
+                // Set unmapped values to zeroGradient. This is the default
+                // action for unmapped fvPatchFields. Note that this bypasses
+                // any special logic for handling unmapped fvPatchFields but
+                // since this is only used inside fvMeshDistribute ...
+
+                tmp<Field<Type>> tfld(bf[patchi].patchInternalField());
+                const Field<Type>& fld = tfld();
+
+                Field<Type> value(bf[patchi]);
+                forAll(directAddressing, i)
+                {
+                    if (directAddressing[i] == -1)
+                    {
+                        value[i] = fld[i];
+                    }
+                }
+                bf[patchi].fvPatchField<Type>::operator=(value);
+            }
+        }
+    }
+
+    return tresult;
+}
+
+
+template<class Type>
+bool Foam::functionObjects::cellDecomposer::mapFieldType() const
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
+
+    const fvMesh& mapRegion =
+        this->mesh_.time().lookupObject<fvMesh>(mapRegion_);
+
+    const labelList patchMap(identity(mapRegion.boundaryMesh().size()));
+
+    const wordList fieldNames
+    (
+        this->mesh_.sortedNames<VolFieldType>(fieldNames_)
+    );
+
+    const bool processed = !fieldNames.empty();
+
+    for (const word& fieldName : fieldNames)
+    {
+        const VolFieldType& field = lookupObject<VolFieldType>(fieldName);
+
+        auto* mapFieldPtr = mapRegion.getObjectPtr<VolFieldType>(fieldName);
+
+        if (!mapFieldPtr)
+        {
+            mapFieldPtr = new VolFieldType
+            (
+                IOobject
+                (
+                    fieldName,
+                    time_.timeName(),
+                    mapRegion,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    IOobject::REGISTER
+                ),
+                mapRegion,
+                dimensioned<Type>(field.dimensions(), Zero)
+            );
+
+            mapFieldPtr->store();
+        }
+
+        auto& mappedField = *mapFieldPtr;
+
+        mappedField = interpolate
+        (
+            field,
+            mapRegion,
+            patchMap,
+            mapPtr_().cellMap(),
+            mapPtr_().faceMap(),
+            false             //allowUnmapped
+        );
+        Log << "    " << fieldName << ": interpolated";
+
+        //evaluateConstraintTypes(mappedField);
+    }
+
+    return processed;
+}
+
+
+template<class Type>
+bool Foam::functionObjects::cellDecomposer::writeFieldType() const
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
+
+    const fvMesh& mapRegion =
+        this->mesh_.time().lookupObject<fvMesh>(mapRegion_);
+
+    const wordList fieldNames
+    (
+        this->mesh_.sortedNames<VolFieldType>(fieldNames_)
+    );
+
+    const bool processed = !fieldNames.empty();
+
+    for (const word& fieldName : fieldNames)
+    {
+        const VolFieldType& mappedField =
+            mapRegion.template lookupObject<VolFieldType>(fieldName);
+
+        mappedField.write();
+
+        Log << "    " << fieldName << ": written";
+    }
+
+    return processed;
+}
+
+
+// ************************************************************************* //

src/functionObjects/field/ddt2/ddt2.C

@@ -172,10 +172,7 @@ bool Foam::functionObjects::ddt2::execute()
 {
     results_.clear();
 
-    wordHashSet candidates
-    (
-        mesh_.names(selectFields_)
-    );
+    wordHashSet candidates(mesh_.names<regIOobject>(selectFields_));
 
     DynamicList<word> missing(selectFields_.size());
     DynamicList<word> ignored(selectFields_.size());

src/functionObjects/field/zeroGradient/zeroGradient.C

@@ -134,10 +134,7 @@ bool Foam::functionObjects::zeroGradient::execute()
 {
     results_.clear();
 
-    wordHashSet candidates
-    (
-        mesh_.names(selectFields_)
-    );
+    wordHashSet candidates(mesh_.names<regIOobject>(selectFields_));
 
     DynamicList<word> missing(selectFields_.size());
     DynamicList<word> ignored(selectFields_.size());

src/functionObjects/forces/propellerInfo/propellerInfo.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021-2022 OpenCFD Ltd.
+    Copyright (C) 2021-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -808,6 +808,7 @@ Foam::functionObjects::propellerInfo::propellerInfo
 )
 :
     forces(name, runTime, dict, false),
+    dict_(dict),
     radius_(0),
     URefPtr_(nullptr),
     rotationMode_(rotationMode::SPECIFIED),
@@ -827,7 +828,8 @@ Foam::functionObjects::propellerInfo::propellerInfo
     interpolationScheme_("cell"),
     wakeFilePtr_(nullptr),
     axialWakeFilePtr_(nullptr),
-    nanValue_(pTraits<scalar>::min)
+    nanValue_(pTraits<scalar>::min),
+    initialised_(false)
 {
     if (readFields)
     {
@@ -855,21 +857,18 @@ bool Foam::functionObjects::propellerInfo::read(const dictionary& dict)
 {
     if (forces::read(dict))
     {
+        dict_ = dict;
+
         radius_ = dict.getScalar("radius");
         URefPtr_.reset(Function1<scalar>::New("URef", dict, &mesh_));
         rotationMode_ = rotationModeNames_.get("rotationMode", dict);
 
-        // Must be set before setting the surface
-        setCoordinateSystem(dict);
-
         writePropellerPerformance_ =
             dict.get<bool>("writePropellerPerformance");
 
         writeWakeFields_ = dict.get<bool>("writeWakeFields");
         if (writeWakeFields_)
         {
-            setSampleDiskSurface(dict);
-
             dict.readIfPresent("interpolationScheme", interpolationScheme_);
 
             dict.readIfPresent("nanValue", nanValue_);
@@ -884,6 +883,19 @@ bool Foam::functionObjects::propellerInfo::read(const dictionary& dict)
 
 bool Foam::functionObjects::propellerInfo::execute()
 {
+    if (!initialised_)
+    {
+        // Must be set before setting the surface
+        setCoordinateSystem(dict_);
+
+        if (writeWakeFields_)
+        {
+            setSampleDiskSurface(dict_);
+        }
+
+        initialised_ = true;
+    }
+
     calcForcesMoments();
 
     createFiles();

src/functionObjects/forces/propellerInfo/propellerInfo.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2021-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -205,6 +205,9 @@ protected:
 
     // Protected data
 
+        //- Copy of dictionary used during construction
+        dictionary dict_;
+
         //- Propeller radius
         scalar radius_;
 
@@ -271,6 +274,9 @@ protected:
             //- scalar::min
             scalar nanValue_;
 
+            //- Initialised flag
+            bool initialised_;
+
 
     // Protected Member Functions
 

src/functionObjects/utilities/Make/files

@@ -39,6 +39,8 @@ systemCall/systemCall.C
 
 timeActivatedFileUpdate/timeActivatedFileUpdate.C
 
+viewFactorHeatFlux/viewFactorHeatFlux.C
+
 writeDictionary/writeDictionary.C
 
 writeObjects/writeObjects.C

src/functionObjects/utilities/viewFactorHeatFlux/viewFactorHeatFlux.C

@@ -0,0 +1,372 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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.
+
+    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 "viewFactorHeatFlux.H"
+#include "addToRunTimeSelectionTable.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+    defineTypeNameAndDebug(viewFactorHeatFlux, 0);
+    addToRunTimeSelectionTable(functionObject, viewFactorHeatFlux, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::functionObjects::viewFactorHeatFlux::initialise()
+{
+    const auto& pbm = mesh_.boundaryMesh();
+    const labelList selectedPatches = pbm.indices("viewFactorWall");
+
+
+    // Initialise output file
+
+    auto& os = file();
+
+    writeHeader(os, "Radiation heat flux");
+    writeCommented(os, "Time");
+
+    for (const label patchi : selectedPatches)
+    {
+        writeTabbed(os, pbm[patchi].name());
+    }
+
+    os  << endl;
+
+
+    // Create compact patch addressing
+
+    label nFace = 0;
+    forAll(selectedPatches, i)
+    {
+        const label patchi = selectedPatches[i];
+        auto slice = compactPatchID_.slice(nFace, pbm[patchi].size());
+        slice = i;
+        nFace += slice.size();
+    }
+
+
+    if (Pstream::parRun())
+    {
+        mapDistPtr_.reset
+        (
+            new IOmapDistribute
+            (
+                IOobject
+                (
+                    "mapDist",
+                    mesh_.facesInstance(),
+                    mesh_,
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE,
+                    IOobject::NO_REGISTER
+                )
+            )
+        );
+
+        auto& mapDist = mapDistPtr_();
+
+        mapDist.distribute(compactPatchID_);
+
+
+        // Convert global addressing into compact form
+
+        labelList compactGlobalIds(mapDist.constructSize(), Zero);
+
+        globalIndex gi(nFace);
+
+        SubList<label>(compactGlobalIds, nFace) = identity(gi.range());
+
+        mapDist.distribute(compactGlobalIds);
+
+        const label nTotalFaces = returnReduce(nFace, sumOp<label>());
+
+        const labelList globalToCompact(invert(nTotalFaces, compactGlobalIds));
+
+        for (labelList& visibleFaces : faceFaces_)
+        {
+            for (label& globalFacei : visibleFaces)
+            {
+                globalFacei = globalToCompact[globalFacei];
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::functionObjects::viewFactorHeatFlux::viewFactorHeatFlux
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict,
+    const bool readFields
+)
+:
+    fvMeshFunctionObject(name, runTime, dict),
+    writeFile(mesh_, name, typeName, dict),
+    qrName_(dict.getOrDefault<word>("qr", "qr")),
+    mapDistPtr_(nullptr),
+    faceFaces_
+    (
+        IOobject
+        (
+            "globalFaceFaces",
+            mesh_.facesInstance(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            IOobject::NO_REGISTER
+        )
+    ),
+    Fij_
+    (
+        IOobject
+        (
+            "F",
+            mesh_.facesInstance(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            IOobject::NO_REGISTER
+        )
+    ),
+    compactPatchID_(faceFaces_.size(), Zero)
+{
+    initialise();
+}
+
+
+Foam::functionObjects::viewFactorHeatFlux::viewFactorHeatFlux
+(
+    const word& name,
+    const objectRegistry& obr,
+    const dictionary& dict,
+    const bool readFields
+)
+:
+    fvMeshFunctionObject(name, obr, dict),
+    writeFile(mesh_, name, typeName, dict),
+    qrName_(dict.getOrDefault<word>("qr", "qr")),
+    mapDistPtr_(nullptr),
+    faceFaces_
+    (
+        IOobject
+        (
+            "globalFaceFaces",
+            mesh_.facesInstance(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            IOobject::NO_REGISTER
+        )
+    ),
+    Fij_
+    (
+        IOobject
+        (
+            "F",
+            mesh_.facesInstance(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            IOobject::NO_REGISTER
+        )
+    ),
+    compactPatchID_(faceFaces_.size(), Zero)
+{
+   initialise();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::functionObjects::viewFactorHeatFlux::read(const dictionary& dict)
+{
+    if (fvMeshFunctionObject::read(dict))
+    {
+        dict.readIfPresent("qr", qrName_);
+
+        return true;
+    }
+
+    return false;
+}
+
+
+bool Foam::functionObjects::viewFactorHeatFlux::execute()
+{
+    return true;
+}
+
+
+bool Foam::functionObjects::viewFactorHeatFlux::write()
+{
+    Log << type() << " " << name() <<  " write:" << nl;
+
+    // Retrieve radiative heat flux field from the mesh database
+
+    const auto* qrPtr = mesh_.cfindObject<volScalarField>(qrName_);
+
+    if (qrPtr)
+    {
+        const auto& qr = *qrPtr;
+
+        // const labelList selectedPatches =
+        //     mesh_.boundaryMesh().indices
+        //     (
+        //         radiation::viewFactor::viewFactorWalls
+        //     );
+        const labelList selectedPatches =
+            mesh_.boundaryMesh().indices("viewFactorWall");
+
+        const auto& pbm = mesh_.boundaryMesh();
+        const label nPatch = selectedPatches.size();
+
+        // Accumulate qr per patch, and
+        // assemble qr contributions in compact addressing
+        scalarList qrPatch(nPatch, Zero);
+        scalarList compactQr(faceFaces_.size());
+
+        label compacti = 0;
+        forAll(selectedPatches, i)
+        {
+            const label patchi = selectedPatches[i];
+            const auto& qrp = qr.boundaryField()[patchi];
+            forAll(qrp, facei)
+            {
+                compactQr[compacti] = qrp[facei];
+                qrPatch[i] += qrp[facei];
+                ++compacti;
+            }
+        }
+
+        reduce(qrPatch, sumOp<scalarList>());
+
+
+        if (Pstream::parRun())
+        {
+            mapDistPtr_->distribute(compactQr);
+        }
+
+        // Populate view factor radiation heat flux matrix
+        scalarSquareMatrix qrVf(nPatch, Zero);
+
+        forAll(Fij_, startFacei)
+        {
+            const scalar qr = compactQr[startFacei];
+            const labelList& visibleSlots = faceFaces_[startFacei];
+            const label i = compactPatchID_[startFacei];
+
+            forAll(visibleSlots, visSloti)
+            {
+                const label sloti = visibleSlots[visSloti];
+                const label j = compactPatchID_[sloti];
+
+                qrVf[i][j] += Fij_[startFacei][visSloti]*qr;
+            }
+        }
+
+        reduce(qrVf, sumOp<scalarSquareMatrix>());
+
+        if (Pstream::master())
+        {
+            Log << "     Writing patch totals to " << file().name()
+                << nl;
+
+            // Write patch sums as time history
+            writeCurrentTime(file());
+
+            for (const auto& qrp : qrPatch)
+            {
+                file() << tab << qrp;
+            }
+
+            file() << endl;
+
+
+            // Write view factor breakdown as matrix per write time
+            auto osPtr = newFileAtTime("viewFactorQr", mesh_.time().value());
+            auto& os = osPtr();
+
+            Log << "     Writing view factor breakdown to " << os.name()
+                << nl;
+
+
+            forAll(selectedPatches, i)
+            {
+                if (i == 0)
+                {
+                    for (const label patchj : selectedPatches)
+                    {
+                        os  << tab << pbm[patchj].name();
+                    }
+                    os  << nl;
+                }
+
+                const label patchi = selectedPatches[i];
+
+                os  << pbm[patchi].name();;
+
+                forAll(selectedPatches, j)
+                {
+                    os  << tab << qrVf[i][j];
+                }
+
+                os  << nl;
+            }
+
+            os  << endl;
+        }
+    }
+
+    Log << endl;
+
+    return true;
+}
+
+
+void Foam::functionObjects::viewFactorHeatFlux::updateMesh(const mapPolyMesh&)
+{
+    NotImplemented;
+}
+
+
+void Foam::functionObjects::viewFactorHeatFlux::movePoints(const polyMesh&)
+{
+    NotImplemented;
+}
+
+
+// ************************************************************************* //

src/functionObjects/utilities/viewFactorHeatFlux/viewFactorHeatFlux.H

@@ -0,0 +1,181 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::functionObjects::viewFactorHeatFlux
+
+Description
+    Determines radiation heat flux between patches when using the viewFactor
+    radiation model.
+
+Usage
+    Minimal example by using \c system/controlDict.functions:
+    \verbatim
+    viewFactorHeatFlux1
+    {
+        // Mandatory entries (unmodifiable)
+        type        viewFactorHeatFlux;
+        libs        (utilityFunctionObjects);
+
+        // Optional entries (runtime modifiable)
+        qr          qr;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property     | Description                        | Type | Req'd | Dflt
+      type         | Type name: viewFactorHeatFlux      | word |  yes  | -
+      libs         | Library name: fieldFunctionObjects | word |  yes  | -
+      qr           | Name of radiation heat flux field  | word |  no   | qr
+    \endtable
+
+See also
+- Foam::radiation::viewFactor
+
+SourceFiles
+    fvMeshFunctionObject.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_functionObjects_viewFactorHeatFlux
+#define Foam_functionObjects_viewFactorHeatFlux
+
+#include "fvMeshFunctionObject.H"
+#include "writeFile.H"
+#include "IOmapDistribute.H"
+#include "labelListIOList.H"
+#include "scalarListIOList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+namespace functionObjects
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class viewFactorHeatFlux Declaration
+\*---------------------------------------------------------------------------*/
+
+class viewFactorHeatFlux
+:
+    public fvMeshFunctionObject,
+    public writeFile
+{
+    // Private Data
+
+        // Name of radiation heat flux field; default = qr
+        word qrName_;
+
+
+        // Fields generated by viewFactorsGen utility
+
+        //- Map
+        autoPtr<IOmapDistribute> mapDistPtr_;
+
+        //- Face addressing
+        labelListIOList faceFaces_;
+
+        //- View factors
+        scalarListIOList Fij_;
+
+        //- Patch indicies in compact addressing
+        scalarList compactPatchID_;
+
+
+    // Private Member Functions
+
+        //- Initialise
+        void initialise();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("viewFactorHeatFlux");
+
+
+    // Constructors
+
+        //- Construct from Time and dictionary
+        viewFactorHeatFlux
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict,
+            const bool readFields = true
+        );
+
+        //- Construct from objectRegistry and dictionary
+        viewFactorHeatFlux
+        (
+            const word& name,
+            const objectRegistry& obr,
+            const dictionary& dict,
+            const bool readFields = true
+        );
+
+        //- No copy construct
+        viewFactorHeatFlux(const viewFactorHeatFlux&) = delete;
+
+        //- No copy assignment
+        void operator=(const viewFactorHeatFlux&) = delete;
+
+
+    //- Destructor
+    virtual ~viewFactorHeatFlux() = default;
+
+
+    // Member Functions
+
+        //- Read the viewFactorHeatFlux data
+        virtual bool read(const dictionary&);
+
+        //- No-op - all operations performed during the call to write()
+        virtual bool execute();
+
+        //- Write the viewFactorHeatFlux per patch
+        virtual bool write();
+
+        //- Update for changes of mesh
+        virtual void updateMesh(const mapPolyMesh&);
+
+        //- Update for mesh point-motion
+        virtual void movePoints(const polyMesh&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/fvMotionSolver/motionDiffusivity/inverseDistance/inverseDistanceDiffusivity.C

@@ -28,7 +28,6 @@ License
 
 #include "inverseDistanceDiffusivity.H"
 #include "addToRunTimeSelectionTable.H"
-#include "patchWave.H"
 #include "HashSet.H"
 #include "surfaceInterpolate.H"
 #include "zeroGradientFvPatchFields.H"

src/fvMotionSolver/motionDiffusivity/inverseVolume/inverseVolumeDiffusivity.C

@@ -27,7 +27,6 @@ License
 
 #include "inverseVolumeDiffusivity.H"
 #include "addToRunTimeSelectionTable.H"
-#include "patchWave.H"
 #include "HashSet.H"
 #include "surfaceInterpolate.H"
 #include "zeroGradientFvPatchFields.H"

src/fvOptions/cellSetOption/cellSetOption.C

@@ -408,9 +408,13 @@ bool Foam::fv::cellSetOption::isActive()
             // Report new volume (if changed)
             setVol();
         }
-        else if (selectionMode_ == smMovingPoints)
+        else if
+        (
+            selectionMode_ == smMovingPoints
+         || selectionMode_ == smGeometric
+        )
         {
-            // Update the cell selection if it moves
+            // Update the cell selection even if mesh is not moving
             setCellSelection();
             setVol();
         }