WIP: INT: Robustness improvements for particle tracking (issue #2537 )

- adjustments similar to openfoam.org changes: * additional protections to prevent division by subnormal numbers and associated floating point errors.
2022-07-19 12:32:43 +02:00
795 changed files with 7249 additions and 19057 deletions
--- a/META-INFO/api-info
+++ b/META-INFO/api-info
@ -1,2 +1,2 @@
 api=2206
-patch=220907
+patch=0
--- a/applications/solvers/acoustic/acousticFoam/Make/options
+++ b/applications/solvers/acoustic/acousticFoam/Make/options
@ -11,5 +11,4 @@ EXE_LIBS = \
    -lfvOptions \
    -lmeshTools \
    -lsampling \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels
--- a/applications/solvers/combustion/fireFoam/Make/options
+++ b/applications/solvers/combustion/fireFoam/Make/options
@ -1,8 +1,8 @@
 EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/finiteArea/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
+    -I${LIB_SRC}/sampling/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
--- a/applications/solvers/compressible/rhoPimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoPimpleFoam/Make/options
@ -26,5 +26,4 @@ EXE_LIBS = \
    -ltopoChangerFvMesh \
    -lsampling \
    -latmosphericModels \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels
--- a/applications/solvers/compressible/rhoSimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoSimpleFoam/Make/options
@ -21,5 +21,4 @@ EXE_LIBS = \
    -lcompressibleTurbulenceModels \
    -lthermoTools \
    -latmosphericModels \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options
@ -28,5 +28,4 @@ EXE_LIBS = \
    -ldynamicMesh \
    -ldynamicFvMesh \
    -ltopoChangerFvMesh \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/Make/options
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/Make/options
@ -22,5 +22,4 @@ EXE_LIBS = \
    -lcompressibleTurbulenceModels \
    -lthermoTools \
    -latmosphericModels \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options
@ -40,5 +40,4 @@ EXE_LIBS = \
    -lthermoTools \
    -lradiationModels \
    -lregionModels \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/Make/options
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/Make/options
@ -36,5 +36,5 @@ EXE_LIBS = \
    -lfaOptions \
    -lregionModels \
    -lsampling \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels
+
--- a/applications/solvers/incompressible/pimpleFoam/Make/options
+++ b/applications/solvers/incompressible/pimpleFoam/Make/options
@ -22,5 +22,4 @@ EXE_LIBS = \
    -ldynamicFvMesh \
    -ltopoChangerFvMesh \
    -latmosphericModels \
-    -lregionFaModels \
-    -lfiniteArea
+    -lregionFaModels
--- a/applications/solvers/lagrangian/coalChemistryFoam/Make/options
+++ b/applications/solvers/lagrangian/coalChemistryFoam/Make/options
@ -1,7 +1,7 @@
 EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/finiteArea/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
--- a/applications/solvers/lagrangian/reactingParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/reactingParcelFoam/Make/options
@ -2,10 +2,10 @@ EXE_INC = \
    -I../reactingParcelFoam \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/finiteArea/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
    -I$(LIB_SRC)/dynamicMesh/lnInclude \
    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
+    -I${LIB_SRC}/sampling/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
--- a/applications/solvers/lagrangian/reactingParcelFoam/reactingHeterogenousParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/reactingParcelFoam/reactingHeterogenousParcelFoam/Make/options
@ -2,8 +2,8 @@ EXE_INC = \
    -I.. \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/finiteArea/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I${LIB_SRC}/sampling/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
    -I$(LIB_SRC)/dynamicMesh/lnInclude \
    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
--- a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/options
@ -1,7 +1,7 @@
 EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/finiteArea/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
    -I$(LIB_SRC)/sampling/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
--- a/applications/solvers/lagrangian/simpleCoalParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/simpleCoalParcelFoam/Make/options
@ -2,7 +2,7 @@ EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/finiteArea/lnInclude \
    -I$(LIB_SRC)/fvOptions/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
--- a/applications/solvers/lagrangian/sprayFoam/Make/options
+++ b/applications/solvers/lagrangian/sprayFoam/Make/options
@ -2,8 +2,8 @@ EXE_INC = \
    -I../reactingParcelFoam \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/finiteArea/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
+    -I${LIB_SRC}/sampling/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
--- a/applications/solvers/lagrangian/sprayFoam/engineFoam/Make/options
+++ b/applications/solvers/lagrangian/sprayFoam/engineFoam/Make/options
@ -4,8 +4,8 @@ EXE_INC = \
    -I../../reactingParcelFoam \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/finiteArea/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
+    -I${LIB_SRC}/sampling/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
--- a/applications/solvers/lagrangian/sprayFoam/simpleSprayFoam/Make/options
+++ b/applications/solvers/lagrangian/sprayFoam/simpleSprayFoam/Make/options
@ -2,7 +2,7 @@ EXE_INC = \
    -I$(FOAM_SOLVERS)/lagrangian/reactingParcelFoam/simpleReactingParcelFoam \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/finiteArea/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
    -I$(LIB_SRC)/sampling/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
--- a/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/Make/options
+++ b/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/Make/options
@ -4,8 +4,8 @@ EXE_INC = \
    -I../../../compressible/rhoPimpleFoam \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/finiteArea/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
+    -I${LIB_SRC}/sampling/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
--- a/applications/solvers/multiphase/MPPICInterFoam/Make/options
+++ b/applications/solvers/multiphase/MPPICInterFoam/Make/options
@ -42,5 +42,4 @@ EXE_LIBS = \
    -lregionModels \
    -lregionFaModels \
    -lsurfaceFilmModels \
-    -lfiniteArea \
    -lcompressibleTwoPhaseMixtureTurbulenceModels
--- a/applications/test/DynamicList2/Test-DynamicList2.C
+++ b/applications/test/DynamicList2/Test-DynamicList2.C
@ -5,7 +5,7 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021-2022 OpenCFD Ltd.
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -104,7 +104,7 @@ void readList
 {
    OTstream os;
    os << input;
-    ITstream is(os.tokens());
+    ITstream is("input", os.tokens());

    is >> output;
 }
@ -118,7 +118,7 @@ void readList
 {
    OTstream os;
    os << input;
-    ITstream is(os.tokens());
+    ITstream is("input", os.tokens());

    is >> output;
 }
--- a/applications/test/FieldFields1/Make/files
+++ b/applications/test/FieldFields1/Make/files
@ -1,3 +0,0 @@
-Test-FieldFields1.C
-
-EXE = $(FOAM_USER_APPBIN)/Test-FieldFields1
--- a/applications/test/FieldFields1/Test-FieldFields1.C
+++ b/applications/test/FieldFields1/Test-FieldFields1.C
@ -1,117 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
-------------------------------------------------------------------------------
-    Copyright (C) 2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Application
-    Test-FieldFields1
-
-\*---------------------------------------------------------------------------*/
-
-#include "symmTensorField.H"
-#include "tensorField.H"
-#include "FieldFields.H"
-#include "Random.H"
-
-using namespace Foam;
-
-
-template<class Cmpt>
-void printFieldField(const FieldField<Field, Cmpt>& ff)
-{
-    forAll(ff, i)
-    {
-        Info<< i << ": " << flatOutput(ff[i]) << nl;
-    }
-    Info<< nl;
-}
-
-
-template<class Type>
-tmp<Field<Type>> randomField(Random& rnd, label dim)
-{
-    auto tfld = tmp<Field<Type>>::New(dim);
-    auto& fld = tfld.ref();
-
-    for (Type& val : fld)
-    {
-        for (direction cmpt=0; cmpt < pTraits<Type>::nComponents; ++cmpt)
-        {
-            setComponent(val, cmpt) = rnd.position<label>(0, 100);
-        }
-    }
-
-    return tfld;
-}
-
-
-template<class Type>
-tmp<Field<Type>> randomField(label dim)
-{
-    Random rnd;
-    return randomField<Type>(rnd, dim);
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    // scalarField
-    {
-        Info<< nl << "scalarFieldField" << nl;
-
-        Random rnd;
-
-        FieldField<Field, scalar> sff1(6);
-        forAll(sff1, i)
-        {
-            sff1.set(i, randomField<scalar>(rnd, 8));
-        }
-
-        printFieldField(sff1);
-
-        Info<< nl << "indexing:" << nl;
-
-        {
-            labelPair index;
-            const label range1 = sff1.size()-1;
-            const label range2 = sff1[0].size()-1;
-
-            for (label iter = 0; iter < 10; ++iter)
-            {
-                index.first() = rnd.position<label>(0, range1);
-                index.second() = rnd.position<label>(0, range2);
-
-                Info<< index << " => " << sff1[index] << nl;
-            }
-        }
-    }
-
-    Info<< nl << "End\n" << nl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //
--- a/applications/test/FieldFields2/Make/files
+++ b/applications/test/FieldFields2/Make/files
@ -1,3 +0,0 @@
-Test-FieldFields2.C
-
-EXE = $(FOAM_USER_APPBIN)/Test-FieldFields2
--- a/applications/test/FieldFields2/Make/options
+++ b/applications/test/FieldFields2/Make/options
@ -1,2 +0,0 @@
-/* EXE_INC = */
-/* EXE_LIBS = */
--- a/applications/test/HashPtrTable/Test-HashPtrTable.C
+++ b/applications/test/HashPtrTable/Test-HashPtrTable.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2017-2022 OpenCFD Ltd.
+    Copyright (C) 2017-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -31,10 +31,7 @@ Description
 #include <memory>
 #include <iostream>
 #include "autoPtr.H"
-#include "refPtr.H"
-#include "tmp.H"
-#include "PtrMap.H"
-#include "primitiveFields.H"
+#include "HashPtrTable.H"

 using namespace Foam;

@ -253,42 +250,6 @@ int main()
        Info<< "Table: " << tbl << nl;
    }

-    {
-        PtrMap<scalarField> fields;
-
-        {
-            const label patchi = 2;
-
-            scalarField fld1(patchi, 5.0);
-            scalarField fld2(patchi, 8.0);
-
-            // assign from tmp<>
-            fields.set( patchi, (fld1 * fld2));
-        }
-
-        {
-            const label patchi = 3;
-
-            scalarField fld1(patchi, 6.0);
-
-            // From tmp (clone)
-            fields.set(patchi, tmp<scalarField>(fld1));
-        }
-
-        {
-            const label patchi = 4;
-
-            // From refPtr
-            fields.set(patchi, refPtr<scalarField>::New(patchi, 10.0));
-        }
-
-        Info<< nl
-            << "PtrMap:" << nl
-            << fields << endl;
-    }
-
-    Info<< "\nEnd" << nl;
-
    return 0;
 }

--- a/applications/test/Hashing1/Make/options
+++ b/applications/test/Hashing1/Make/options
@ -1,3 +1,3 @@
-EXE_INC = $(c++LESSWARN)
+EXE_INC = ${c++LESSWARN}

 /* EXE_LIBS = */
--- a/applications/test/Hashing2/Make/options
+++ b/applications/test/Hashing2/Make/options
@ -1,3 +1,3 @@
-EXE_INC = $(c++LESSWARN)
+EXE_INC = ${c++LESSWARN}

 /* EXE_LIBS = */
--- a/applications/test/HashingSpeed/Make/options
+++ b/applications/test/HashingSpeed/Make/options
@ -1,3 +1,3 @@
-EXE_INC = $(c++LESSWARN)
+EXE_INC = ${c++LESSWARN}

 /* EXE_LIBS = */
--- a/applications/test/Polynomial/Test-Polynomial.C
+++ b/applications/test/Polynomial/Test-Polynomial.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2020-2022 OpenCFD Ltd.
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -132,7 +132,7 @@ int main(int argc, char *argv[])
        OTstream os;
        os << poly;

-        ITstream is(std::move(os.tokens()));
+        ITstream is("input", std::move(os.tokens()));
        is >> polyfunc;
    }

--- a/applications/test/TestTools/TestTools.H
+++ b/applications/test/TestTools/TestTools.H
@ -28,14 +28,15 @@ Description

 \*---------------------------------------------------------------------------*/

-#include "scalar.H"
+using namespace Foam;
+
+#include "floatScalar.H"
+#include "doubleScalar.H"
 #include "complex.H"
 #include "Matrix.H"
 #include "Random.H"
 #include <chrono>

-using namespace Foam;
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

 // Total number of unit tests
--- a/applications/test/autoPtr/Test-autoPtr.C
+++ b/applications/test/autoPtr/Test-autoPtr.C
@ -5,7 +5,7 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -27,6 +27,7 @@ License

 // #define Foam_autoPtr_deprecate_setMethod

+#include <memory>
 #include "autoPtr.H"
 #include "labelList.H"
 #include "ListOps.H"
@ -76,14 +77,6 @@ autoPtr<T> testNullReturn2()
 }


-template<class T>
-struct DerivedList : public List<T>
-{
-    // Inherit constructors
-    using List<T>::List;
-};
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 //  Main program:

@ -261,16 +254,6 @@ int main(int argc, char *argv[])
        // autoPtr<labelList> ptr2 = testNullReturn2<labelList>();
    }

-    {
-        auto input1 = autoPtr<DerivedList<label>>::New(label(10), 1);
-        auto input2 = autoPtr<DerivedList<scalar>>::New(label(10), 1.0);
-
-        autoPtr<labelList> ptr1(std::move(input1));
-
-        // Does not compile: ptr1 = std::move(input2);
-        // Does not compile: ptr1 = autoPtr<List<scalar>>::New(label(10), 2);
-    }
-

    return 0;
 }
--- a/applications/test/dimField/Test-dimField.C
+++ b/applications/test/dimField/Test-dimField.C
@ -42,7 +42,7 @@ namespace Foam
    // Something like an internal state field. Probably only dimensionless
    typedef DimensionedField<uint8_t, volMesh> dimUint8Field;

-    defineTemplateTypeNameAndDebug(dimUint8Field, 0);
+    defineTemplate2TypeNameAndDebug(dimUint8Field, 0);

 } // End namespace Foam
 #endif
--- a/applications/test/dynamicLibrary/Test-dynamicLibrary.C
+++ b/applications/test/dynamicLibrary/Test-dynamicLibrary.C
@ -5,7 +5,7 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020-2022 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -46,8 +46,8 @@ int main(int argc, char *argv[])
    profiling::disable(); // No profiling output
    argList::noBanner();
    argList::noParallel();
-    argList::noFunctionObjects();
    argList::removeOption("case");
+    argList::removeOption("noFunctionObjects");
    argList::addBoolOption("no-close", "Skip dlclose");
    argList::addBoolOption("quiet", "Disable verbosity");

--- a/applications/test/faceHashing/Make/options
+++ b/applications/test/faceHashing/Make/options
@ -1,3 +1,3 @@
-EXE_INC = $(c++LESSWARN)
+EXE_INC = ${c++LESSWARN}

 /* EXE_LIBS = */
--- a/applications/test/fieldTypes/Make/options
+++ b/applications/test/fieldTypes/Make/options
@ -1,7 +1,7 @@
 EXE_INC = \
    -I$(LIB_SRC)/finiteArea/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I${LIB_SRC}/finiteVolume/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \

 EXE_LIBS = \
    -lfiniteArea \
--- a/applications/test/fvc2D/Test-fvc2D.C
+++ b/applications/test/fvc2D/Test-fvc2D.C
@ -39,7 +39,7 @@ namespace Foam
 {
    typedef GeometricField<vector2D, fvPatchField, volMesh> volVector2DField;

-    defineTemplateTypeNameAndDebug(volVector2DField::Internal, 0);
+    defineTemplate2TypeNameAndDebug(volVector2DField::Internal, 0);
    defineTemplateTypeNameAndDebug(volVector2DField, 0);

    typedef fvPatchField<vector2D> fvPatchVector2DField;
--- a/applications/test/matrices/DiagonalMatrix/Test-DiagonalMatrix.C
+++ b/applications/test/matrices/DiagonalMatrix/Test-DiagonalMatrix.C
@ -38,7 +38,8 @@ Description

 #include "DiagonalMatrix.H"
 #include "RectangularMatrix.H"
-#include "scalar.H"
+#include "floatScalar.H"
+#include "doubleScalar.H"
 #include "complex.H"
 #include "TestTools.H"

--- a/applications/test/matrices/RectangularMatrix/Test-RectangularMatrix.C
+++ b/applications/test/matrices/RectangularMatrix/Test-RectangularMatrix.C
@ -44,7 +44,8 @@ Note

 #include "RectangularMatrix.H"
 #include "SquareMatrix.H"
-#include "scalar.H"
+#include "floatScalar.H"
+#include "doubleScalar.H"
 #include "complex.H"
 #include "IOmanip.H"
 #include "TestTools.H"
--- a/applications/test/matrices/SquareMatrix/Test-SquareMatrix.C
+++ b/applications/test/matrices/SquareMatrix/Test-SquareMatrix.C
@ -45,7 +45,8 @@ Note
 #include "scalarMatrices.H"
 #include "RectangularMatrix.H"
 #include "SquareMatrix.H"
-#include "scalar.H"
+#include "floatScalar.H"
+#include "doubleScalar.H"
 #include "complex.H"
 #include "IOmanip.H"
 #include "TestTools.H"
--- a/applications/test/matrices/SymmetricSquareMatrix/Test-SymmetricSquareMatrix.C
+++ b/applications/test/matrices/SymmetricSquareMatrix/Test-SymmetricSquareMatrix.C
@ -46,7 +46,8 @@ Note
 #include "RectangularMatrix.H"
 #include "SquareMatrix.H"
 #include "SymmetricSquareMatrix.H"
-#include "scalar.H"
+#include "floatScalar.H"
+#include "doubleScalar.H"
 #include "complex.H"
 #include "IOmanip.H"
 #include "Random.H"
--- a/applications/test/openmp/Make/options
+++ b/applications/test/openmp/Make/options
@ -1,4 +1,4 @@
-EXE_INC = $(COMP_OPENMP) /* -UUSE_OMP */
+EXE_INC = ${COMP_OPENMP} /* -UUSE_OMP */

 /* Mostly do not need to explicitly link openmp libraries */
-/* EXE_LIBS = $(LINK_OPENMP) */
+/* EXE_LIBS = ${LINK_OPENMP} */
--- a/applications/test/parallel-chunks/Test-parallel-chunks.C
+++ b/applications/test/parallel-chunks/Test-parallel-chunks.C
@ -55,7 +55,7 @@ void do_exchangeBuf
    const bool wait
 )
 {
-    const label startOfRequests = UPstream::nRequests();
+    const label startOfRequests = Pstream::nRequests();

    // Set up receives
    // ~~~~~~~~~~~~~~~
@ -132,7 +132,7 @@ void do_exchangeContainer
    const bool wait
 )
 {
-    const label startOfRequests = UPstream::nRequests();
+    const label startOfRequests = Pstream::nRequests();

    // Set up receives
    // ~~~~~~~~~~~~~~~
--- a/applications/test/refPtr/Test-refPtr.C
+++ b/applications/test/refPtr/Test-refPtr.C
@ -5,7 +5,7 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020-2022 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@ -19,9 +19,6 @@ Description
 \*---------------------------------------------------------------------------*/

 #include "primitiveFields.H"
-#include "autoPtr.H"
-#include "refPtr.H"
-#include "tmp.H"
 #include "Switch.H"

 using namespace Foam;
@ -32,36 +29,6 @@ struct myScalarField : public scalarField
 };


-template<class T>
-void constructInfo()
-{
-    Info<< " move-constructible:"
-        << std::is_move_constructible<T>::value
-        << " move-assignable:"
-        << std::is_move_assignable<T>::value
-        << " nothrow:"
-        << std::is_nothrow_move_assignable<T>::value
-        << " trivially:"
-        << std::is_trivially_move_assignable<T>::value
-        << nl;
-}
-
-
-template<class T>
-void printInfo(const autoPtr<T>& item, const bool verbose = false)
-{
-    Info<< "autoPtr good:" << Switch::name(item.good())
-        << " addr: " << Foam::name(item.get());
-
-    constructInfo<autoPtr<T>>();
-
-    if (verbose && item)
-    {
-        Info<< "content: " << item() << nl;
-    }
-}
-
-
 template<class T>
 void printInfo(const refPtr<T>& item, const bool verbose = false)
 {
@ -70,24 +37,15 @@ void printInfo(const refPtr<T>& item, const bool verbose = false)
        << " addr: " << Foam::name(item.get())
        << " movable:" << Switch(item.movable());

-    constructInfo<refPtr<T>>();
-
-    if (verbose && item)
-    {
-        Info<< "content: " << item() << nl;
-    }
-}
-
-
-template<class T>
-void printInfo(const tmp<T>& item, const bool verbose = false)
-{
-    Info<< "tmp good:" << Switch::name(item.good())
-        << " pointer:" << Switch::name(item.is_pointer())
-        << " addr: " << Foam::name(item.get())
-        << " movable:" << Switch(item.movable());
-
-    constructInfo<tmp<T>>();
+    Info<< " move-constructible:"
+        << std::is_move_constructible<refPtr<T>>::value
+        << " move-assignable:"
+        << std::is_move_assignable<refPtr<T>>::value
+        << " nothrow:"
+        << std::is_nothrow_move_assignable<refPtr<T>>::value
+        << " trivially:"
+        << std::is_trivially_move_assignable<refPtr<T>>::value
+        << nl;

    if (verbose && item)
    {
@ -143,62 +101,6 @@ int main()
        printInfo(tfld3, true);
    }

-    {
-        refPtr<scalarField> tfld1;
-        auto aptr = autoPtr<scalarField>::New(2, scalar(2));
-
-        tmp<scalarField> tfld2;
-        printInfo(tfld2, true);
-
-        tfld2 = new scalarField(10, Zero);
-
-    /*
-        tfld2 = aptr.get();
-
-        // tfld1.reset(aptr);
-        // tfld1 = std::move(aptr);
-        // tfld1 = aptr;
-
-        Info<< nl << "From autoPtr" << nl;
-        printInfo(aptr, true);
-        //& Info<< nl << "Construct from autoPtr" << nl;
-        //& // refPtr<scalarField> tfld2(autoPtr<scalarField>::New(10, scalar(2)));
-        //& printInfo(tfld2, true);
-*/
-    }
-
-    {
-        auto aptr1 = autoPtr<labelField>::New(2, Zero);
-        //auto aptr1 = autoPtr<scalarField>::New(2, scalar(2));
-        auto aptr2 = autoPtr<scalarField>::New(2, scalar(2));
-
-        refPtr<scalarField> tfld2(std::move(aptr2));
-
-        // aptr2 = std::move(aptr1);
-    }
-
-    {
-        auto tptr1 = tmp<labelField>::New(2, Zero);
-        auto aptr1 = autoPtr<labelField>::New(2, Zero);
-        auto tfld2 = refPtr<labelField>::New(2, Zero);
-
-        // Deleted: refPtr<labelField> tfld1(aptr1);
-        refPtr<labelField> tfld1;
-
-        // refPtr<labelField> tfld1(std::move(tptr1));
-        // refPtr<labelField> tfld1(tptr1);
-
-        tfld1 = std::move(aptr1);
-
-        // tfld1.reset(aptr1);
-        // tfld1.reset(tfld2);
-
-        // tfld1 = std::move(tptr1);
-        // Does not compile: tfld1.ref(tptr1);
-        // Deleted: tfld1.cref(tptr1);
-        // Deleted: tfld1.ref(aptr1);
-    }
-
    Info<< "\nEnd" << endl;

    return 0;
--- a/applications/test/scalarOps/Test-scalarOps.C
+++ b/applications/test/scalarOps/Test-scalarOps.C
@ -126,9 +126,6 @@ int main(int argc, char *argv[])
    testDivide<vector>(vectors);


-    //(void) compareOp<vector>()(vector::zero, vector::one);
-
-
    Info<< "\nEnd\n" << endl;

    return 0;
--- a/applications/test/tensorFieldFields1/Make/files
+++ b/applications/test/tensorFieldFields1/Make/files
@ -0,0 +1,3 @@
+Test-tensorFieldFields1.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-tensorFieldFields1
--- a/applications/test/tensorFieldFields1/Make/options
+++ b/applications/test/tensorFieldFields1/Make/options
--- a/applications/test/tensorFieldFields1/Test-tensorFieldFields1.C
+++ b/applications/test/tensorFieldFields1/Test-tensorFieldFields1.C
@ -24,7 +24,7 @@ License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

 Application
-    Test-FieldFields2
+    Test-tensorFieldFields1

 \*---------------------------------------------------------------------------*/

@ -73,29 +73,21 @@ void allocComponents
 }


-template<class Type>
-tmp<Field<Type>> randomField(Random& rnd, label dim)
+vectorField randomVectorField(label size)
 {
-    auto tfld = tmp<Field<Type>>::New(dim);
-    auto& fld = tfld.ref();
+    Random rnd;

-    for (Type& val : fld)
+    vectorField vf(size);
+
+    forAll(vf, i)
    {
-        for (direction cmpt=0; cmpt < pTraits<Type>::nComponents; ++cmpt)
+        for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt)
        {
-            setComponent(val, cmpt) = rnd.position<label>(0, 100);
+            vf[i][cmpt] = rnd.position<label>(0, 100);
        }
    }

-    return tfld;
-}
-
-
-template<class Type>
-tmp<Field<Type>> randomField(label dim)
-{
-    Random rnd;
-    return randomField<Type>(rnd, dim);
+    return vf;
 }


@ -199,7 +191,7 @@ int main(int argc, char *argv[])
        printFieldField(sf1);

        Info<< nl;
-        for (direction cmpt = 0; cmpt < pTraits<vector>::nComponents; ++cmpt)
+        for (direction cmpt = 0; cmpt < vector::nComponents; ++cmpt)
        {
            unzipRow(sf1, cmpt, slice[0]);

@ -261,9 +253,8 @@ int main(int argc, char *argv[])
    {
        Info<< nl << "vectorField" << nl;

-        Random rnd;
        FieldField<Field, vector> vf1(1);
-        vf1.set(0, randomField<vector>(rnd, 4));
+        vf1.set(0, new vectorField(randomVectorField(4)));

        FixedList<FieldField<Field, scalar>, 3> cmpts;
        allocComponents(cmpts, 4);
--- a/applications/test/tensorFields1/Test-tensorFields1.C
+++ b/applications/test/tensorFields1/Test-tensorFields1.C
@ -34,29 +34,21 @@ Application

 using namespace Foam;

-template<class Type>
-tmp<Field<Type>> randomField(Random& rnd, label dim)
+vectorField randomVectorField(label size)
 {
-    auto tfld = tmp<Field<Type>>::New(dim);
-    auto& fld = tfld.ref();
+    Random rnd;

-    for (Type& val : fld)
+    vectorField vf(size);
+
+    forAll(vf, i)
    {
-        for (direction cmpt=0; cmpt < pTraits<Type>::nComponents; ++cmpt)
+        for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt)
        {
-            setComponent(val, cmpt) = rnd.position<label>(0, 100);
+            vf[i][cmpt] = rnd.position<label>(0, 100);
        }
    }

-    return tfld;
-}
-
-
-template<class Type>
-tmp<Field<Type>> randomField(label dim)
-{
-    Random rnd;
-    return randomField<Type>(rnd, dim);
+    return vf;
 }


@ -259,7 +251,7 @@ int main(int argc, char *argv[])
    {
        Info<< nl << "vectorField" << nl;

-        vectorField vf1(randomField<vector>(4));
+        vectorField vf1(randomVectorField(4));
        FixedList<scalarField, 3> cmpts(scalarField(vf1.size()));

        Info<< nl
--- a/applications/test/timeSelector/Test-timeSelector.C
+++ b/applications/test/timeSelector/Test-timeSelector.C
@ -63,6 +63,7 @@ int main(int argc, char *argv[])
    argList::addNote("Test timeSelector and TimePaths");

    timeSelector::addOptions(true, true);
+    argList::noLibs();
    argList::noFunctionObjects();

    argList::addOption("relative", "PATH", "Test relativePath");
--- a/applications/test/tmp/Test-tmp.C
+++ b/applications/test/tmp/Test-tmp.C
@ -19,9 +19,6 @@ Description
 \*---------------------------------------------------------------------------*/

 #include "primitiveFields.H"
-#include "autoPtr.H"
-#include "refPtr.H"
-#include "tmp.H"
 #include "Switch.H"

 using namespace Foam;
@ -129,28 +126,6 @@ int main()
        printInfo(tfld2);
    }

-    {
-        auto tptr1 = refPtr<labelField>::New(2, Zero);
-        auto aptr1 = autoPtr<labelField>::New(2, Zero);
-
-        // Deleted: tmp<labelField> tfld1(aptr1);
-        // tmp<labelField> tfld1(std::move(aptr1));
-        // tmp<labelField> tfld1(std::move(tptr1));
-        tmp<labelField> tfld1;
-        //tfld1.cref(tptr1);
-        //tfld1.cref(aptr1);
-
-        // refPtr<labelField> tfld1(std::move(tptr1));
-        // refPtr<labelField> tfld1(tptr1);
-
-        // tfld1 = std::move(aptr1);
-
-        // tfld1 = std::move(tptr1);
-        // Does not compile: tfld1.ref(tptr1);
-        // Deleted: tfld1.cref(tptr1);
-        // Deleted: tfld1.ref(aptr1);
-    }
-
    Info<< "\nEnd" << endl;

    return 0;
--- a/applications/test/volField/Test-volField.C
+++ b/applications/test/volField/Test-volField.C
@ -42,7 +42,7 @@ namespace Foam
    // - still needs some basic boundary conditions!!
    typedef GeometricField<uint8_t, fvPatchField, volMesh> volUint8Field;

-    defineTemplateTypeNameAndDebug(volUint8Field, 0);
+    defineTemplate2TypeNameAndDebug(volUint8Field, 0);

 } // End namespace Foam
 #endif
--- a/applications/utilities/finiteArea/checkFaMesh/printMeshSummary.H
+++ b/applications/utilities/finiteArea/checkFaMesh/printMeshSummary.H
@ -116,16 +116,12 @@ Description
    {
        const faPatch& p = patches[patchi];

-        // Report physical size (nEdges) not virtual size
-        Info<< "  " << "patch " << p.index()
-            << " (size: " << returnReduce(p.nEdges(), sumOp<label>())
-            << ") name: " << p.name()
-            << nl;
+         Info<< "  " << "patch " << p.index()
+             << " (size: " << returnReduce(p.size(), sumOp<label>())
+             << ") name: " << p.name()
+             << nl;
    }

-    Info<< "----------------" << nl
-        << "Used polyPatches: " << flatOutput(aMesh.whichPolyPatches()) << nl;
-

    // Geometry information
    Info<< nl;
--- a/applications/utilities/finiteArea/makeFaMesh/printMeshSummary.H
+++ b/applications/utilities/finiteArea/makeFaMesh/printMeshSummary.H
@ -5,7 +5,7 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021-2022 OpenCFD Ltd.
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@ -36,17 +36,12 @@ Description
    {
        const faPatch& p = patches[patchi];

-        // Report physical size (nEdges) not virtual size
-        Info<< "  " << "patch " << p.index()
-            << " (size: " << returnReduce(p.nEdges(), sumOp<label>())
-            << ") name: " << p.name()
-            << nl;
+         Info<< "  " << "patch " << p.index()
+             << " (size: " << returnReduce(p.size(), sumOp<label>())
+             << ") name: " << p.name()
+             << nl;
    }

-    Info<< "----------------" << nl
-        << "Used polyPatches: " << flatOutput(aMesh.whichPolyPatches()) << nl;
-
-
    // Geometry information
    Info<< nl;
    {
--- a/applications/utilities/mesh/advanced/refinementLevel/refinementLevel.C
+++ b/applications/utilities/mesh/advanced/refinementLevel/refinementLevel.C
@ -226,18 +226,14 @@ int main(int argc, char *argv[])
    // Add the boundary patches
    const polyBoundaryMesh& patches = mesh.boundaryMesh();

-    polyPatchList newPatches(patches.size());
+    List<polyPatch*> p(patches.size());

-    forAll(newPatches, patchi)
+    forAll(p, patchi)
    {
-        newPatches.set
-        (
-            patchi,
-            patches[patchi].clone(fMesh.boundaryMesh())
-        );
+        p[patchi] = patches[patchi].clone(fMesh.boundaryMesh()).ptr();
    }

-    fMesh.addFvPatches(newPatches);
+    fMesh.addFvPatches(p);


    // Refinement level
--- a/applications/utilities/mesh/advanced/snappyRefineMesh/snappyRefineMesh.C
+++ b/applications/utilities/mesh/advanced/snappyRefineMesh/snappyRefineMesh.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018-2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -205,43 +205,34 @@ label addPatch(polyMesh& mesh, const word& patchName)
    {
        const polyBoundaryMesh& patches = mesh.boundaryMesh();

-        polyPatchList newPatches(patches.size() + 1);
-
-        label nPatches = 0;
+        List<polyPatch*> newPatches(patches.size() + 1);

        // Add empty patch as 0th entry (Note: only since subsetMesh wants this)
        patchi = 0;
-        {
-            newPatches.set
-            (
-                nPatches,
-                new emptyPolyPatch
-                (
-                    patchName,
-                    0,
-                    mesh.nInternalFaces(),
-                    nPatches,
-                    patches,
-                    emptyPolyPatch::typeName
-                )
-            );
-            ++nPatches;
-        }

-        for (const polyPatch& pp : patches)
-        {
-            newPatches.set
+        newPatches[patchi] =
+            new emptyPolyPatch
            (
-                nPatches,
+                Foam::word(patchName),
+                0,
+                mesh.nInternalFaces(),
+                patchi,
+                patches,
+                emptyPolyPatch::typeName
+            );
+
+        forAll(patches, i)
+        {
+            const polyPatch& pp = patches[i];
+
+            newPatches[i+1] =
                pp.clone
                (
                    patches,
-                    nPatches,
+                    i+1,
                    pp.size(),
                    pp.start()
-                )
-            );
-            ++nPatches;
+                ).ptr();
        }

        mesh.removeBoundary();
--- a/applications/utilities/mesh/conversion/gmshToFoam/gmshToFoam.C
+++ b/applications/utilities/mesh/conversion/gmshToFoam/gmshToFoam.C
@ -1678,7 +1678,7 @@ int main(int argc, char *argv[])
    label defaultPatchID = mesh.boundaryMesh().findPatchID(defaultFacesName);
    if (mesh.boundaryMesh()[defaultPatchID].size() == 0)
    {
-        polyPatchList newPatches((mesh.boundaryMesh().size() - 1));
+        List<polyPatch*> newPatchPtrList((mesh.boundaryMesh().size() - 1));
        label newPatchi = 0;
        forAll(mesh.boundaryMesh(), patchi)
        {
@ -1686,21 +1686,18 @@ int main(int argc, char *argv[])
            {
                const polyPatch& patch = mesh.boundaryMesh()[patchi];

-                newPatches.set
+                newPatchPtrList[newPatchi] = patch.clone
                (
+                    mesh.boundaryMesh(),
                    newPatchi,
-                    patch.clone
-                    (
-                        mesh.boundaryMesh(),
-                        newPatchi,
-                        patch.size(),
-                        patch.start()
-                    )
-                );
-                ++newPatchi;
+                    patch.size(),
+                    patch.start()
+                ).ptr();
+
+                newPatchi++;
            }
        }
-        repatcher.changePatches(newPatches);
+        repatcher.changePatches(newPatchPtrList);
    }

    // Set the precision of the points data to 10
--- a/applications/utilities/mesh/generation/blockMesh/mergePatchPairs.H
+++ b/applications/utilities/mesh/generation/blockMesh/mergePatchPairs.H
@ -1,56 +1,30 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
-------------------------------------------------------------------------------
-    Copyright (C) 2009 OpenFOAM Foundation
-    Copyright (C) 2017-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
-
-Description
-    Handle merging of patch pairs
-
-\*---------------------------------------------------------------------------*/
-
+// Handle merging of patch pairs
 {
    wordPairList mergePatchPairs;

-    // Read in a list of merge patch pairs
+    // Read in a list of dictionaries for the merge patch pairs
    if
    (
        meshDict.readIfPresent("mergePatchPairs", mergePatchPairs)
     && mergePatchPairs.size()
    )
    {
-        Info<< "Merging " << mergePatchPairs.size() << " patch pairs" << nl;
+        Info<< "Creating merge patch pairs" << nl << endl;

-        // Cleanup
-        wordHashSet cleanupPatches(4*mergePatchPairs.size());
-        wordHashSet cleanupPointZones(2*mergePatchPairs.size());
-        wordHashSet cleanupFaceZones(2*mergePatchPairs.size());
-
-        Info<< "    Adding point and face zones" << endl;
+        Info<< "Adding point and face zones" << endl;
        {
-            const auto& pbm = mesh.boundaryMesh();
-
-            auto& pzs = mesh.pointZones(); pzs.clearAddressing();
-            auto& fzs = mesh.faceZones();  fzs.clearAddressing();
+            auto& pzs = mesh.pointZones();
+            pzs.clearAddressing();
+            auto& fzs = mesh.faceZones();
+            fzs.clearAddressing();

            forAll(mergePatchPairs, pairi)
            {
-                // Patch pairs
-                const polyPatch& patch0 = pbm[mergePatchPairs[pairi].first()];
-                const polyPatch& patch1 = pbm[mergePatchPairs[pairi].second()];
-
                const word mergeName
                (
                    mergePatchPairs[pairi].first()
                  + mergePatchPairs[pairi].second()
-                  + Foam::name(pairi)
+                  + name(pairi)
                );

                // An empty zone for cut points
@ -63,35 +37,40 @@ Description
                        pzs
                    )
                );
-                cleanupPointZones.insert(pzs.last().name());

-                // Coupling side 0 (master)
+                // Master patch
+                const word masterPatchName(mergePatchPairs[pairi].first());
+                const polyPatch& masterPatch =
+                    mesh.boundaryMesh()[masterPatchName];
+
                fzs.append
                (
                    new faceZone
                    (
-                        mergeName + "Side0Zone",
-                        identity(patch0.range()),
+                        mergeName + "MasterZone",
+                        identity(masterPatch.range()),
                        false, // none are flipped
                        fzs.size(),
                        fzs
                    )
                );
-                cleanupFaceZones.insert(fzs.last().name());

-                // Coupling side 1 (slave)
+                // Slave patch
+                const word slavePatchName(mergePatchPairs[pairi].second());
+                const polyPatch& slavePatch =
+                    mesh.boundaryMesh()[slavePatchName];
+
                fzs.append
                (
                    new faceZone
                    (
-                        mergeName + "Side1Zone",
-                        identity(patch1.range()),
+                        mergeName + "SlaveZone",
+                        identity(slavePatch.range()),
                        false, // none are flipped
                        fzs.size(),
                        fzs
                    )
                );
-                cleanupFaceZones.insert(fzs.last().name());

                // An empty zone for cut faces
                fzs.append
@ -103,38 +82,35 @@ Description
                        fzs
                    )
                );
-                cleanupFaceZones.insert(fzs.last().name());
-            }
+            }  // end of all merge pairs
        }


-        Info<< "    Merging with attachPolyTopoChanger" << endl;
+
+        Info<< "Creating attachPolyTopoChanger" << endl;
        attachPolyTopoChanger polyMeshAttacher(mesh);
-        polyMeshAttacher.resize(1);
+        polyMeshAttacher.setSize(mergePatchPairs.size());

        forAll(mergePatchPairs, pairi)
        {
-            cleanupPatches.insert(mergePatchPairs[pairi].first());
-            cleanupPatches.insert(mergePatchPairs[pairi].second());
-
            const word mergeName
            (
                mergePatchPairs[pairi].first()
              + mergePatchPairs[pairi].second()
-              + Foam::name(pairi)
+              + name(pairi)
            );

            // Add the sliding interface mesh modifier
            polyMeshAttacher.set
            (
-                0,
+                pairi,
                new slidingInterface
                (
-                    "couple" + Foam::name(pairi),
+                    "couple" + name(pairi),
                    pairi,
                    polyMeshAttacher,
-                    mergeName + "Side0Zone",
-                    mergeName + "Side1Zone",
+                    mergeName + "MasterZone",
+                    mergeName + "SlaveZone",
                    mergeName + "CutPointZone",
                    mergeName + "CutFaceZone",
                    mergePatchPairs[pairi].first(),
@ -144,135 +120,12 @@ Description
                    intersection::VISIBLE
                )
            );
-
-            polyMeshAttacher.attach(false);  // Do not yet remove empty patches
        }

-        // Re-do the boundary patches, removing empty merge patches
-        // but keeping any other empty patches
-        {
-            const polyBoundaryMesh& oldPatches = mesh.boundaryMesh();
-
-            polyPatchList newPatches(oldPatches.size());
-            label nNewPatches = 0;
-
-            wordHashSet removedPatches(cleanupPatches.capacity());
-
-            forAll(oldPatches, patchi)
-            {
-                const word& patchName = oldPatches[patchi].name();
-
-                if
-                (
-                    !cleanupPatches.found(patchName)
-                 || returnReduceOr(oldPatches[patchi].size())
-                )
-                {
-                    newPatches.set
-                    (
-                        nNewPatches,
-                        oldPatches[patchi].clone
-                        (
-                            mesh.boundaryMesh(),
-                            nNewPatches,
-                            oldPatches[patchi].size(),
-                            oldPatches[patchi].start()
-                        )
-                    );
-
-                    ++nNewPatches;
-                }
-                else
-                {
-                    removedPatches.insert(patchName);
-                }
-            }
-
-            newPatches.resize(nNewPatches);
-
-            mesh.removeBoundary();
-            mesh.addPatches(newPatches);
-
-            Info<< "Removed " << removedPatches.size()
-                << " empty merged patches:" << nl
-                << "    " << flatOutput(removedPatches.sortedToc()) << endl;
-        }
-
-        // Cleanup empty merged point zones
-        {
-            PtrList<pointZone>& zones = mesh.pointZones();
-            mesh.pointZones().clearAddressing();
-
-            wordHashSet removedZones(2*zones.size());
-
-            label nZones = 0;
-            forAll(zones, zonei)
-            {
-                if
-                (
-                    !cleanupPointZones.found(zones[zonei].name())
-                 || returnReduceOr(zones[zonei].size())
-                )
-                {
-                    zones.set(nZones, zones.release(zonei));
-                    zones[nZones].index() = nZones;  // re-index
-                    ++nZones;
-                }
-                else
-                {
-                    removedZones.insert(zones[zonei].name());
-                }
-            }
-            zones.resize(nZones);
-
-            if (removedZones.size())
-            {
-                Info<< "Removed " << removedZones.size()
-                    << " empty point zones:" << nl
-                    << "    " << flatOutput(removedZones.sortedToc()) << endl;
-            }
-        }
-
-        // Cleanup empty merged face zones
-        {
-            PtrList<faceZone>& zones = mesh.faceZones();
-            mesh.faceZones().clearAddressing();
-
-            wordHashSet removedZones(2*zones.size());
-
-            label nZones = 0;
-            forAll(zones, zonei)
-            {
-                if
-                (
-                    !cleanupFaceZones.found(zones[zonei].name())
-                 || returnReduceOr(zones[zonei].size())
-                )
-                {
-                    zones.set(nZones, zones.release(zonei));
-                    zones[nZones].index() = nZones;  // re-index
-                    ++nZones;
-                }
-                else
-                {
-                    removedZones.insert(zones[zonei].name());
-                }
-            }
-            zones.resize(nZones);
-
-            if (removedZones.size())
-            {
-                Info<< "Removed " << removedZones.size()
-                    << " empty merged face zones:" << nl
-                    << "    " << flatOutput(removedZones.sortedToc()) << endl;
-            }
-        }
+        polyMeshAttacher.attach(true);
    }
    else
    {
-        Info<< "No patch pairs to merge" << endl;
+        Info<< nl << "There are no merge patch pairs" << endl;
    }
 }
-
-
-// ************************************************************************* //
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/PrintTable/PrintTable.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/PrintTable/PrintTable.C
@ -152,7 +152,7 @@ void Foam::PrintTable<KeyType, DataType>::print

        os  << endl;

-        const List<KeyType> sortedTable(combinedTable.sortedToc());
+        const List<KeyType>& sortedTable = combinedTable.sortedToc();

        forAll(sortedTable, keyI)
        {
--- a/applications/utilities/mesh/generation/foamyMesh/foamyQuadMesh/foamyQuadMesh.C
+++ b/applications/utilities/mesh/generation/foamyMesh/foamyQuadMesh/foamyQuadMesh.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2013-2016 OpenFOAM Foundation
-    Copyright (C) 2021-2022 OpenCFD Ltd.
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -167,32 +167,28 @@ int main(int argc, char *argv[])
    );

    Info<< "Constructing patches." << endl;
-    polyPatchList newPatches(poly2DMesh.patchNames().size());
-    label nPatches = 0;
+    List<polyPatch*> patches(poly2DMesh.patchNames().size());
+    label countPatches = 0;

-    forAll(newPatches, patchi)
+    forAll(patches, patchi)
    {
        if (poly2DMesh.patchSizes()[patchi] != 0)
        {
-            newPatches.set
+            patches[countPatches] = new polyPatch
            (
-                nPatches,
-                new polyPatch
-                (
-                    poly2DMesh.patchNames()[patchi],
-                    poly2DMesh.patchSizes()[patchi],
-                    poly2DMesh.patchStarts()[patchi],
-                    nPatches,
-                    pMesh.boundaryMesh(),
-                    word::null
-                )
+                poly2DMesh.patchNames()[patchi],
+                poly2DMesh.patchSizes()[patchi],
+                poly2DMesh.patchStarts()[patchi],
+                countPatches,
+                pMesh.boundaryMesh(),
+                word::null
            );

-            ++nPatches;
+            countPatches++;
        }
    }
-    newPatches.resize(nPatches);
-    pMesh.addPatches(newPatches);
+    patches.setSize(countPatches);
+    pMesh.addPatches(patches);

    if (extrude)
    {
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
@ -1342,17 +1342,6 @@ int main(int argc, char *argv[])
    const snapParameters snapParams(snapDict, dryRun);


-    Info<< "Setting refinement level of surface to be consistent"
-        << " with curvature." << endl;
-    surfaces.setCurvatureMinLevelFields
-    (
-        refineParams.curvature(),
-        meshRefiner.meshCutter().level0EdgeLength()
-    );
-    Info<< "Checked curvature refinement in = "
-        << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
-
-

    // Add all the cellZones and faceZones
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
--- a/applications/utilities/mesh/manipulation/autoPatch/autoPatch.C
+++ b/applications/utilities/mesh/manipulation/autoPatch/autoPatch.C
@ -191,7 +191,7 @@ int main(int argc, char *argv[])
    const PtrList<boundaryPatch>& patches = bMesh.patches();

    // Create new list of patches with old ones first
-    polyPatchList newPatches(patches.size());
+    List<polyPatch*> newPatchPtrList(patches.size());

    newPatchi = 0;

@ -200,41 +200,34 @@ int main(int argc, char *argv[])
    {
        const polyPatch& patch = mesh.boundaryMesh()[patchi];

-        newPatches.set
-        (
-            newPatchi,
+        newPatchPtrList[newPatchi] =
            patch.clone
            (
                mesh.boundaryMesh(),
                newPatchi,
                patch.size(),
                patch.start()
-            )
-        );
+            ).ptr();

-        ++newPatchi;
+        newPatchi++;
    }

    // Add new ones with empty size.
    for (label patchi = newPatchi; patchi < patches.size(); patchi++)
    {
-        const word& patchName = patches[patchi].name();
+        const boundaryPatch& bp = patches[patchi];

-        newPatches.set
+        newPatchPtrList[newPatchi] = polyPatch::New
        (
+            polyPatch::typeName,
+            bp.name(),
+            0,
+            mesh.nFaces(),
            newPatchi,
-            polyPatch::New
-            (
-                polyPatch::typeName,
-                patchName,
-                0,
-                mesh.nFaces(),
-                newPatchi,
-                mesh.boundaryMesh()
-            )
-        );
+            mesh.boundaryMesh()
+        ).ptr();

-        ++newPatchi;
+        newPatchi++;
    }

    if (!overwrite)
@ -245,7 +238,7 @@ int main(int argc, char *argv[])

    // Change patches
    repatchPolyTopoChanger polyMeshRepatcher(mesh);
-    polyMeshRepatcher.changePatches(newPatches);
+    polyMeshRepatcher.changePatches(newPatchPtrList);


    // Change face ordering
--- a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
@ -40,7 +40,6 @@ License
 #include "vtkSurfaceWriter.H"
 #include "checkTools.H"
 #include "treeBoundBox.H"
-#include "syncTools.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

@ -49,113 +48,25 @@ void Foam::checkPatch
 (
    const bool allGeometry,
    const word& name,
-    const polyMesh& mesh,
    const PatchType& pp,
-    const labelList& meshFaces,
-    const labelList& meshEdges,
    pointSet& points
 )
 {
-    typedef typename PatchType::surfaceTopo TopoType;
-
-    const label globalSize = returnReduce(pp.size(), sumOp<label>());
-
    Info<< "    "
        << setw(20) << name
-        << setw(9) << globalSize
+        << setw(9) << returnReduce(pp.size(), sumOp<label>())
        << setw(9) << returnReduce(pp.nPoints(), sumOp<label>());

-    if (globalSize == 0)
-    {
-        Info<< setw(34) << "ok (empty)";
-    }
-    else if (Pstream::parRun())
-    {
-        // Parallel - use mesh edges
-        // - no check for point-pinch
-        // - no check for consistent orientation (if that is posible to
-        //   check?)
-
-        // (see addPatchCellLayer::globalEdgeFaces)
-        // From mesh edge to global face labels. Non-empty sublists only for
-        // pp edges.
-        labelListList globalEdgeFaces(mesh.nEdges());
-
-        const labelListList& edgeFaces = pp.edgeFaces();
-
-        // Global numbering
-        const globalIndex globalFaces(mesh.nFaces());
-
-        forAll(edgeFaces, edgei)
-        {
-            label meshEdgei = meshEdges[edgei];
-            const labelList& eFaces = edgeFaces[edgei];
-
-            // Store face and processor as unique tag.
-            labelList& globalEFaces = globalEdgeFaces[meshEdgei];
-            globalEFaces.setSize(eFaces.size());
-            forAll(eFaces, i)
-            {
-                globalEFaces[i] = globalFaces.toGlobal(meshFaces[eFaces[i]]);
-            }
-            //Pout<< "At edge:" << meshEdgei
-            //    << " ctr:" << mesh.edges()[meshEdgei].centre(mesh.points())
-            //    << " have eFaces:" << globalEdgeFaces[meshEdgei]
-            //    << endl;
-        }
-
-        //DebugVar(globalEdgeFaces);
-
-
-        // Synchronise across coupled edges.
-        syncTools::syncEdgeList
-        (
-            mesh,
-            globalEdgeFaces,
-            ListOps::uniqueEqOp<label>(),
-            labelList()             // null value
-        );
-
-        //DebugVar(globalEdgeFaces);
-        
-        label labelTyp = TopoType::MANIFOLD;
-        forAll(meshEdges, edgei)
-        {
-            const label meshEdgei = meshEdges[edgei];
-            const labelList& globalEFaces = globalEdgeFaces[meshEdgei];
-            if (globalEFaces.size() == 1)
-            {
-                //points.insert(mesh.edges()[meshEdgei]);
-                labelTyp = max(labelTyp, TopoType::OPEN);
-            }
-            else if (globalEFaces.size() == 0 || globalEFaces.size() > 2)
-            {
-                points.insert(mesh.edges()[meshEdgei]);
-                labelTyp = max(labelTyp, TopoType::ILLEGAL);
-            }
-        }
-        reduce(labelTyp, maxOp<label>());
-
-        if (labelTyp == TopoType::MANIFOLD)
-        {
-            Info<< setw(34) << "ok (closed singly connected)";
-        }
-        else if (labelTyp == TopoType::OPEN)
-        {
-            Info<< setw(34)
-                << "ok (non-closed singly connected)";
-        }
-        else
-        {
-            Info<< setw(34)
-                << "multiply connected (shared edge)";
-        }
-    }
-    else
+    if (!Pstream::parRun())
    {
+        typedef typename PatchType::surfaceTopo TopoType;
        TopoType pTyp = pp.surfaceType();

-        if (pTyp == TopoType::MANIFOLD)
+        if (pp.empty())
+        {
+            Info<< setw(34) << "ok (empty)";
+        }
+        else if (pTyp == TopoType::MANIFOLD)
        {
            if (pp.checkPointManifold(true, &points))
            {
@ -638,8 +549,15 @@ Foam::label Foam::checkTopology
    );

    {
-        Info<< "\nChecking patch topology for multiply connected"
+        if (!Pstream::parRun())
+        {
+            Info<< "\nChecking patch topology for multiply connected"
                << " surfaces..." << endl;
+        }
+        else
+        {
+            Info<< "\nChecking basic patch addressing..." << endl;
+        }

        const polyBoundaryMesh& patches = mesh.boundaryMesh();

@ -648,8 +566,11 @@ Foam::label Foam::checkTopology
        Info<< "    "
            << setw(20) << "Patch"
            << setw(9) << "Faces"
-            << setw(9) << "Points"
-            << "Surface topology";
+            << setw(9) << "Points";
+        if (!Pstream::parRun())
+        {
+            Info<< setw(34) << "Surface topology";
+        }
        if (allGeometry)
        {
            Info<< " Bounding box";
@ -662,16 +583,7 @@ Foam::label Foam::checkTopology

            if (!isA<processorPolyPatch>(pp))
            {
-                checkPatch
-                (
-                    allGeometry,
-                    pp.name(),
-                    mesh,
-                    pp,
-                    identity(pp.size(), pp.start()),
-                    pp.meshEdges(),
-                    points
-                );
+                checkPatch(allGeometry, pp.name(), pp, points);
                Info<< endl;
            }
        }
@ -680,8 +592,15 @@ Foam::label Foam::checkTopology
    }

    {
-        Info<< "\nChecking faceZone topology for multiply connected"
-            << " surfaces..." << endl;
+        if (!Pstream::parRun())
+        {
+            Info<< "\nChecking faceZone topology for multiply connected"
+                << " surfaces..." << endl;
+        }
+        else
+        {
+            Info<< "\nChecking basic faceZone addressing..." << endl;
+        }

        Pout.setf(ios_base::left);

@ -692,8 +611,12 @@ Foam::label Foam::checkTopology
            Info<< "    "
                << setw(20) << "FaceZone"
                << setw(9) << "Faces"
-                << setw(9) << "Points"
-                << setw(34) << "Surface topology";
+                << setw(9) << "Points";
+
+            if (!Pstream::parRun())
+            {
+                Info<< setw(34) << "Surface topology";
+            }
            if (allGeometry)
            {
                Info<< " Bounding box";
@ -702,16 +625,7 @@ Foam::label Foam::checkTopology

            for (const faceZone& fz : faceZones)
            {
-                checkPatch
-                (
-                    allGeometry,
-                    fz.name(),
-                    mesh,
-                    fz(),           // patch
-                    fz,             // mesh face labels
-                    fz.meshEdges(), // mesh edge labels
-                    points
-                );
+                checkPatch(allGeometry, fz.name(), fz(), points);
                Info<< endl;
            }
        }
--- a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.H
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.H
@ -1,4 +1,4 @@
-#include "labelList.H"
+#include "label.H"
 #include "autoPtr.H"

 namespace Foam
@ -14,10 +14,7 @@ namespace Foam
    (
        const bool allGeometry,
        const word& name,
-        const polyMesh& mesh,
        const PatchType& pp,
-        const labelList& meshFaces,
-        const labelList& meshEdges,
        pointSet& points
    );

--- a/applications/utilities/mesh/manipulation/mergeMeshes/mergePolyMesh.C
+++ b/applications/utilities/mesh/manipulation/mergeMeshes/mergePolyMesh.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2016-2022 OpenCFD Ltd.
+    Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -136,55 +136,46 @@ void Foam::mergePolyMesh::sortProcessorPatches()

    const polyBoundaryMesh& oldPatches = boundaryMesh();

-    polyPatchList newPatches(oldPatches.size());
+    DynamicList<polyPatch*> newPatches(oldPatches.size());

    labelList oldToSorted(oldPatches.size());

-    label nPatches = 0;
-
    forAll(oldPatches, patchi)
    {
        const polyPatch& pp = oldPatches[patchi];

        if (!isA<processorPolyPatch>(pp))
        {
-            newPatches.set
+            oldToSorted[patchi] = newPatches.size();
+            newPatches.append
            (
-                nPatches,
-                pp.clone
-                (
-                    oldPatches,
-                    nPatches,
-                    0,
-                    nInternalFaces()
-                )
-            );
-
-            oldToSorted[patchi] = nPatches;
-            ++nPatches;
-        }
-    }
-
-    forAll(oldPatches, patchi)
-    {
-        const polyPatch& pp = oldPatches[patchi];
-
-        if (isA<processorPolyPatch>(pp))
-        {
-            newPatches.set
-            (
-                nPatches,
                pp.clone
                (
                    oldPatches,
                    oldToSorted[patchi],
                    0,
                    nInternalFaces()
-                )
+                ).ptr()
            );
+        }
+    }
+    forAll(oldPatches, patchi)
+    {
+        const polyPatch& pp = oldPatches[patchi];

-            oldToSorted[patchi] = nPatches;
-            ++nPatches;
+        if (isA<processorPolyPatch>(pp))
+        {
+            oldToSorted[patchi] = newPatches.size();
+            newPatches.append
+            (
+                pp.clone
+                (
+                    oldPatches,
+                    oldToSorted[patchi],
+                    0,
+                    nInternalFaces()
+                ).ptr()
+            );
        }
    }

@ -470,7 +461,7 @@ void Foam::mergePolyMesh::merge()
    {
        Info<< "Copying old patches" << endl;

-        polyPatchList newPatches(patchNames_.size());
+        List<polyPatch*> newPatches(patchNames_.size());

        const polyBoundaryMesh& oldPatches = boundaryMesh();

@ -479,11 +470,7 @@ void Foam::mergePolyMesh::merge()

        for (patchi = 0; patchi < oldPatches.size(); patchi++)
        {
-            newPatches.set
-            (
-                patchi,
-                oldPatches[patchi].clone(oldPatches)
-            );
+            newPatches[patchi] = oldPatches[patchi].clone(oldPatches).ptr();
        }

        Info<< "Adding new patches. " << endl;
@ -500,18 +487,15 @@ void Foam::mergePolyMesh::merge()
            dict.set("nFaces", 0);
            dict.set("startFace", endOfLastPatch);

-            newPatches.set
+            newPatches[patchi] =
            (
-                patchi,
+                polyPatch::New
                (
-                    polyPatch::New
-                    (
-                        patchNames_[patchi],
-                        dict,
-                        patchi,
-                        oldPatches
-                    )
-                )
+                    patchNames_[patchi],
+                    dict,
+                    patchi,
+                    oldPatches
+                ).ptr()
            );
        }

--- a/applications/utilities/mesh/manipulation/mirrorMesh/mirrorFvMesh.C
+++ b/applications/utilities/mesh/manipulation/mirrorMesh/mirrorFvMesh.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2017-2022 OpenCFD Ltd.
+    Copyright (C) 2017-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -404,24 +404,20 @@ Foam::mirrorFvMesh::mirrorFvMesh
    fvMesh& pMesh = mirrorMeshPtr_();

    // Add the boundary patches
-    polyPatchList newPatches(newPatchSizes.size());
+    List<polyPatch*> p(newPatchSizes.size());

-    forAll(newPatches, patchi)
+    forAll(p, patchi)
    {
-        newPatches.set
+        p[patchi] = boundaryMesh()[patchi].clone
        (
+            pMesh.boundaryMesh(),
            patchi,
-            boundaryMesh()[patchi].clone
-            (
-                pMesh.boundaryMesh(),
-                patchi,
-                newPatchSizes[patchi],
-                newPatchStarts[patchi]
-            )
-        );
+            newPatchSizes[patchi],
+            newPatchStarts[patchi]
+        ).ptr();
    }

-    pMesh.addPatches(newPatches);
+    pMesh.addPatches(p);
 }


--- a/applications/utilities/mesh/manipulation/renumberMesh/Make/options
+++ b/applications/utilities/mesh/manipulation/renumberMesh/Make/options
@ -3,7 +3,7 @@
 */

 EXE_INC = \
-    $(COMP_FLAGS) \
+    ${COMP_FLAGS} \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
    -I$(LIB_SRC)/dynamicMesh/lnInclude \
--- a/applications/utilities/mesh/manipulation/stitchMesh/stitchMesh.C
+++ b/applications/utilities/mesh/manipulation/stitchMesh/stitchMesh.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2017-2022 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -564,7 +564,7 @@ int main(int argc, char *argv[])

            mesh.faceZones()
            (
-                mergePatchName + "Side0Zone",
+                mergePatchName + "MasterZone",
                true // verbose
            ).resetAddressing(std::move(faceIds), false);

@ -574,7 +574,7 @@ int main(int argc, char *argv[])

            mesh.faceZones()
            (
-                mergePatchName + "Side1Zone",
+                mergePatchName + "SlaveZone",
                true // verbose
            ).resetAddressing(std::move(faceIds), false);

@ -595,8 +595,8 @@ int main(int argc, char *argv[])
                    "couple" + Foam::name(actioni),
                    0,
                    stitcher,
-                    mergePatchName + "Side0Zone",
-                    mergePatchName + "Side1Zone",
+                    mergePatchName + "MasterZone",
+                    mergePatchName + "SlaveZone",
                    mergePatchName + "CutPointZone",
                    cutZoneName,
                    masterPatchName,
--- a/applications/utilities/mesh/manipulation/transformPoints/transformPoints.C
+++ b/applications/utilities/mesh/manipulation/transformPoints/transformPoints.C
@ -62,8 +62,8 @@ Usage
    -rotate-z angle
        Rotate (degrees) about z-axis.

-     or -yawPitchRoll : (yaw pitch roll) degrees
-     or -rollPitchYaw : (roll pitch yaw) degrees
+     or -yawPitchRoll (yawdegrees pitchdegrees rolldegrees)
+     or -rollPitchYaw (rolldegrees pitchdegrees yawdegrees)

    -scale scalar|vector
        Scale the points by the given scalar or vector on output.
@ -268,18 +268,15 @@ int main(int argc, char *argv[])
    );
    argList::addBoolOption
    (
-        "auto-centre",
-        "Use bounding box centre as centre for rotations"
+        "auto-origin",
+        "Use bounding box centre as origin for rotations"
    );
    argList::addOption
    (
-        "centre",
+        "origin",
        "point",
-        "Use specified <point> as centre for rotations"
+        "Use specified <point> as origin for rotations"
    );
-    argList::addOptionCompat("auto-centre", {"auto-origin", 2206});
-    argList::addOptionCompat("centre", {"origin", 2206});
-
    argList::addOption
    (
        "rotate",
@ -440,18 +437,18 @@ int main(int argc, char *argv[])
            points += v;
        }

-        vector rotationCentre;
-        bool useRotationCentre = args.readIfPresent("centre", rotationCentre);
-        if (args.found("auto-centre") && !useRotationCentre)
+        vector origin;
+        bool useOrigin = args.readIfPresent("origin", origin);
+        if (args.found("auto-origin") && !useOrigin)
        {
-            useRotationCentre = true;
-            rotationCentre = boundBox(points).centre();
+            useOrigin = true;
+            origin = boundBox(points).centre();
        }

-        if (useRotationCentre)
+        if (useOrigin)
        {
-            Info<< "Set centre of rotation to " << rotationCentre << endl;
-            points -= rotationCentre;
+            Info<< "Set origin for rotations to " << origin << endl;
+            points -= origin;
        }


@ -548,15 +545,15 @@ int main(int argc, char *argv[])
            }
        }

-        if (useRotationCentre)
-        {
-            Info<< "Unset centre of rotation from " << rotationCentre << endl;
-            points += rotationCentre;
-        }
-
        // Output scaling
        applyScaling(points, getScalingOpt("scale", args));

+        if (useOrigin)
+        {
+            Info<< "Unset origin for rotations from " << origin << endl;
+            points += origin;
+        }
+

        // Set the precision of the points data to 10
        IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
--- a/applications/utilities/miscellaneous/OSspecific/addr2line/addr2line.C
+++ b/applications/utilities/miscellaneous/OSspecific/addr2line/addr2line.C
@ -43,7 +43,7 @@ Description
 static void usage();
 static void version();
 static std::string getLine(const std::string&, const std::string&);
-static std::string pipeOpen(const std::string& cmd, const int lineNum = 0);
+static std::string pOpen(const std::string&, int line=0);


 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -132,50 +132,51 @@ void version()
 }


-// Read up to and including lineNum from the piped command
-// Return the final line read
-std::string pipeOpen(const std::string& cmd, int lineNum)
+std::string pOpen(const std::string& cmd, int line)
 {
-    std::string str;
+    std::string res;

-    FILE* handle = popen(cmd.c_str(), "r");
-    if (!handle) return str;
+    FILE* cmdPipe = popen(cmd.c_str(), "r");
+    if (!cmdPipe) return res;

    char* buf = nullptr;
-    size_t len = 0;
-    ssize_t nread;

-    // Read lineNum number of lines
-    for
-    (
-        int cnt = 0;
-        cnt <= lineNum && (nread = ::getline(&buf, &len, handle)) >= 0;
-        ++cnt
-    )
+    // Read line number of lines
+    for (int cnt = 0; cnt <= line; ++cnt)
    {
-        if (cnt == lineNum)
-        {
-            // Retain the last line, trimming trailing newline
-            str.assign(buf);
+        size_t linecap = 0;
+        ssize_t linelen = ::getline(&buf, &linecap, cmdPipe);

-            if (str.size())
+        if (linelen < 0)
+        {
+            break;
+        }
+
+        if (cnt == line)
+        {
+            res = std::string(buf);
+
+            // Trim trailing newline
+            if (res.size())
            {
-                str.resize(str.size()-1);
+                res.resize(res.size()-1);
            }
+            break;
        }
    }

-    free(buf);
-    pclose(handle);
+    if (buf) free(buf);

-    return str;
+    pclose(cmdPipe);
+
+    return res;
 }


 std::string getLine(const std::string& filename, const std::string& addr)
 {
    std::string line =
-        pipeOpen
+        pOpen
        (
            "echo 'image lookup -va " + addr
          + "'"
--- a/applications/utilities/miscellaneous/foamHasLibrary/foamHasLibrary.C
+++ b/applications/utilities/miscellaneous/foamHasLibrary/foamHasLibrary.C
@ -5,7 +5,7 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020-2022 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -65,8 +65,8 @@ int main(int argc, char *argv[])
    profiling::disable(); // No profiling output
    argList::noBanner();
    argList::noParallel();
-    argList::noFunctionObjects();
    argList::removeOption("case");
+    argList::removeOption("noFunctionObjects");
    argList::addBoolOption
    (
        "or",
--- a/applications/utilities/parallelProcessing/decomposePar/decomposePar.C
+++ b/applications/utilities/parallelProcessing/decomposePar/decomposePar.C
@ -1014,9 +1014,7 @@ int main(int argc, char *argv[])
                                Time::controlDictName,
                                args.rootPath(),
                                args.caseName()
-                              / ("processor" + Foam::name(proci)),
-                                args.allowFunctionObjects(),
-                                args.allowLibs()
+                              / ("processor" + Foam::name(proci))
                            )
                        );
                    }
--- a/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C
+++ b/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C
@ -244,9 +244,7 @@ int main(int argc, char *argv[])
            (
                Time::controlDictName,
                args.rootPath(),
-                args.caseName()/("processor" + Foam::name(proci)),
-                args.allowFunctionObjects(),
-                args.allowLibs()
+                args.caseName()/("processor" + Foam::name(proci))
            )
        );
    }
--- a/applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C
+++ b/applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C
@ -731,9 +731,7 @@ int main(int argc, char *argv[])
            (
                Time::controlDictName,
                args.rootPath(),
-                args.caseName()/("processor" + Foam::name(proci)),
-                args.allowFunctionObjects(),
-                args.allowLibs()
+                args.caseName()/("processor" + Foam::name(proci))
            )
        );
    }
@ -1186,7 +1184,7 @@ int main(int argc, char *argv[])
                        false
                    );

-                    masterMeshPtr.cref(fvMeshes[0]);
+                    masterMeshPtr = fvMeshes[0];
                }


--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C
@ -168,6 +168,7 @@ int main(int argc, char *argv[])

    // Less frequently used - reduce some clutter
    argList::setAdvanced("decomposeParDict");
+    argList::setAdvanced("noFunctionObjects");

    argList::addVerboseOption("Additional verbosity");

--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/writeDimFields.H
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/writeDimFields.H
@ -65,7 +65,7 @@ label writeDimFields
 {
    typedef VolumeInternalField<Type> FieldType;

-    const auto& mesh = refCast<const fvMesh>(ensMesh.mesh());
+    const fvMesh& mesh = dynamicCast<const fvMesh>(ensMesh.mesh());

    label count = 0;

--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/writeVolFields.H
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/writeVolFields.H
@ -88,7 +88,7 @@ label writeVolFields
 {
    typedef VolumeField<Type> FieldType;

-    const auto& mesh = refCast<const fvMesh>(ensMesh.mesh());
+    const fvMesh& mesh = dynamicCast<const fvMesh>(ensMesh.mesh());

    label count = 0;

--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK.C
@ -232,13 +232,13 @@ vtk::outputOptions getOutputOptions(const argList& args)

        if (!args.found("ascii"))
        {
-            if (sizeof(float) != 4 || sizeof(label) != 4)
+            if (sizeof(floatScalar) != 4 || sizeof(label) != 4)
            {
                opts.ascii(true);

                WarningInFunction
                    << "Using ASCII rather than legacy binary VTK format since "
-                    << "float and/or label are not 4 bytes in size."
+                    << "floatScalar and/or label are not 4 bytes in size."
                    << nl << endl;
            }
            else
@ -268,6 +268,7 @@ int main(int argc, char *argv[])

    // Less frequently used - reduce some clutter
    argList::setAdvanced("decomposeParDict");
+    argList::setAdvanced("noFunctionObjects");

    argList::addVerboseOption("Additional verbosity");

--- a/applications/utilities/postProcessing/lagrangian/particleTracks/particleTracks.C
+++ b/applications/utilities/postProcessing/lagrangian/particleTracks/particleTracks.C
@ -87,6 +87,7 @@ int main(int argc, char *argv[])

    // Less frequently used - reduce some clutter
    argList::setAdvanced("decomposeParDict");
+    argList::setAdvanced("noFunctionObjects");

    argList::addOption
    (
--- a/applications/utilities/postProcessing/miscellaneous/profilingSummary/profilingSummary.C
+++ b/applications/utilities/postProcessing/miscellaneous/profilingSummary/profilingSummary.C
@ -5,7 +5,7 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2022 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -100,9 +100,7 @@ int main(int argc, char *argv[])
            (
                Time::controlDictName,
                args.rootPath(),
-                args.caseName()/("processor" + Foam::name(proci)),
-                args.allowFunctionObjects(),
-                args.allowLibs()
+                args.caseName()/("processor" + Foam::name(proci))
            )
        );
    }
--- a/applications/utilities/preProcessing/changeDictionary/changeDictionary.C
+++ b/applications/utilities/preProcessing/changeDictionary/changeDictionary.C
@ -662,8 +662,8 @@ int main(int argc, char *argv[])
                // Read dictionary
                // Note: disable class type checking so we can load field
                Info<< "Loading dictionary " << fieldName << endl;
-                const word oldTypeName = localIOdictionary::typeName;
-                const_cast<word&>(localIOdictionary::typeName) = word::null;
+                const word oldTypeName = IOdictionary::typeName;
+                const_cast<word&>(IOdictionary::typeName) = word::null;

                IOobject fieldHeader
                (
@ -675,15 +675,11 @@ int main(int argc, char *argv[])
                    false
                );

-                if (fieldHeader.typeHeaderOk<localIOdictionary>(false))
+                if (fieldHeader.typeHeaderOk<IOdictionary>(false))
                {
-                    //IOdictionary fieldDict(fieldHeader);
-                    //- dictionaries to-be-changed are either boundary
-                    //  or field dictionary. Both are processor-local.
-                    localIOdictionary fieldDict(fieldHeader);
+                    IOdictionary fieldDict(fieldHeader);

-                    const_cast<word&>(localIOdictionary::typeName) =
-                        oldTypeName;
+                    const_cast<word&>(IOdictionary::typeName) = oldTypeName;

                    // Fake type back to what was in field
                    const_cast<word&>(fieldDict.type()) =
--- a/applications/utilities/preProcessing/faceAgglomerate/faceAgglomerate.C
+++ b/applications/utilities/preProcessing/faceAgglomerate/faceAgglomerate.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2020,2022 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -96,15 +96,7 @@ int main(int argc, char *argv[])

    label nCoarseFaces = 0;

-
-    const auto& patchesDict =
-        agglomDict.optionalSubDict
-        (
-            "patchAgglomeration",
-            keyType::LITERAL
-        );
-
-    for (const entry& dEntry : patchesDict)
+    for (const entry& dEntry : agglomDict)
    {
        labelList patchids = boundary.indices(dEntry.keyword());

@ -120,7 +112,7 @@ int main(int argc, char *argv[])
                (
                    pp.localFaces(),
                    pp.localPoints(),
-                    dEntry.dict()
+                    agglomDict.subDict(pp.name())
                );

                agglomObject.agglomerate();
--- a/applications/utilities/preProcessing/mapFieldsPar/MapVolFields.H
+++ b/applications/utilities/preProcessing/mapFieldsPar/MapVolFields.H
@ -43,8 +43,7 @@ void evaluateConstraintTypes(GeometricField<Type, fvPatchField, volMesh>& fld)
 {
    auto& fldBf = fld.boundaryFieldRef();

-    const UPstream::commsTypes commsType = UPstream::defaultCommsType;
-    const label startOfRequests = UPstream::nRequests();
+    const UPstream::commsTypes commsType(UPstream::defaultCommsType);

    if
    (
@ -52,6 +51,8 @@ void evaluateConstraintTypes(GeometricField<Type, fvPatchField, volMesh>& fld)
     || commsType == UPstream::commsTypes::nonBlocking
    )
    {
+        const label startOfRequests = UPstream::nRequests();
+
        forAll(fldBf, patchi)
        {
            fvPatchField<Type>& tgtField = fldBf[patchi];
@ -67,7 +68,11 @@ void evaluateConstraintTypes(GeometricField<Type, fvPatchField, volMesh>& fld)
        }

        // Wait for outstanding requests
-        if (commsType == UPstream::commsTypes::nonBlocking)
+        if
+        (
+            UPstream::parRun()
+         && commsType == UPstream::commsTypes::nonBlocking
+        )
        {
            UPstream::waitRequests(startOfRequests);
        }
--- a/applications/utilities/preProcessing/setExprBoundaryFields/setExprBoundaryFields.C
+++ b/applications/utilities/preProcessing/setExprBoundaryFields/setExprBoundaryFields.C
@ -56,7 +56,6 @@ using namespace Foam;

 int main(int argc, char *argv[])
 {
-    // Normally without functionObjects, with -withFunctionObjects to enable
    argList::noFunctionObjects(true);

    // No -constant, no special treatment for 0/
--- a/applications/utilities/preProcessing/setExprFields/setExprFields.C
+++ b/applications/utilities/preProcessing/setExprFields/setExprFields.C
@ -494,7 +494,6 @@ void evaluate

 int main(int argc, char *argv[])
 {
-    // Normally without functionObjects, with -withFunctionObjects to enable
    argList::noFunctionObjects(true);

    // No -constant, no special treatment for 0/
@ -742,7 +741,7 @@ int main(int argc, char *argv[])
            );
        }

-        if (args.allowFunctionObjects())
+        if (args.found("withFunctionObjects"))
        {
            runTime.functionObjects().start();
        }
--- a/applications/utilities/preProcessing/setFields/Make/options
+++ b/applications/utilities/preProcessing/setFields/Make/options
@ -1,10 +1,8 @@
 EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/finiteArea/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude

 EXE_LIBS = \
    -lfiniteVolume \
-    -lfiniteArea \
    -lmeshTools \
    -lgenericPatchFields
--- a/applications/utilities/preProcessing/setFields/setFields.C
+++ b/applications/utilities/preProcessing/setFields/setFields.C
@ -6,7 +6,6 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -38,312 +37,215 @@ Description
 #include "argList.H"
 #include "Time.H"
 #include "fvMesh.H"
-#include "faMesh.H"
 #include "topoSetSource.H"
 #include "cellSet.H"
 #include "faceSet.H"
 #include "volFields.H"
-#include "areaFields.H"

 using namespace Foam;

-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-//- Simple tuple of field type and name (read from stream)
-class fieldDescription
-{
-    word type_;
-    word name_;
-
-public:
-
-    const word& type() const noexcept { return type_; }
-    const word& name() const noexcept { return name_; }
-
-    explicit fieldDescription(Istream& is)
-    {
-        is >> type_;
-        is >> name_;
-
-        // Eg, read as "volScalarFieldValue", but change to "volScalarField"
-        if (type_.ends_with("Value"))
-        {
-            type_.erase(type_.size()-5);
-        }
-    }
-};
-
-
-// Consume unused field information
-template<class Type>
-bool consumeUnusedType(const fieldDescription& fieldDesc, Istream& is)
-{
-    typedef GeometricField<Type, faPatchField, areaMesh> fieldType1;
-    typedef GeometricField<Type, fvPatchField, volMesh>  fieldType2;
-    //? typedef GeometricField<Type, faePatchField, areaMesh> fieldType3;
-    //? typedef GeometricField<Type, fvsPatchField, volMesh>  fieldType4;
-
-    if
-    (
-        fieldDesc.type() == fieldType1::typeName
-     || fieldDesc.type() == fieldType2::typeName
-    )
-    {
-        (void) pTraits<Type>(is);
-        return true;
-    }
-
-    return false;
-}
-
-
-// Consume unused field information
-static bool consumeUnused(const fieldDescription& fieldDesc, Istream& is)
-{
-    return
-    (
-        consumeUnusedType<scalar>(fieldDesc, is)
-     || consumeUnusedType<vector>(fieldDesc, is)
-     || consumeUnusedType<sphericalTensor>(fieldDesc, is)
-     || consumeUnusedType<symmTensor>(fieldDesc, is)
-     || consumeUnusedType<tensor>(fieldDesc, is)
-    );
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-// Setting volume fields
 template<class Type>
 bool setCellFieldType
 (
-    const fieldDescription& fieldDesc,
+    const word& fieldTypeDesc,
    const fvMesh& mesh,
    const labelList& selectedCells,
-    Istream& is
+    Istream& fieldValueStream
 )
 {
    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;

-    if (fieldDesc.type() != fieldType::typeName)
+    if (fieldTypeDesc != fieldType::typeName + "Value")
    {
        return false;
    }

-    // Get value from stream
-    const Type fieldValue = pTraits<Type>(is);
-
+    word fieldName(fieldValueStream);

    // Check the current time directory
    IOobject fieldHeader
    (
-        fieldDesc.name(),
-        mesh.thisDb().time().timeName(),
-        mesh.thisDb(),
+        fieldName,
+        mesh.time().timeName(),
+        mesh,
        IOobject::MUST_READ
    );

-    bool found = fieldHeader.typeHeaderOk<fieldType>(true);
-
-    if (!found)
+    // Check the "constant" directory
+    if (!fieldHeader.typeHeaderOk<fieldType>(true))
    {
-        // Fallback to "constant" directory
        fieldHeader = IOobject
        (
-            fieldDesc.name(),
-            mesh.thisDb().time().constant(),
-            mesh.thisDb(),
+            fieldName,
+            mesh.time().constant(),
+            mesh,
            IOobject::MUST_READ
        );
-        found = fieldHeader.typeHeaderOk<fieldType>(true);
    }

-    // Field exists
-    if (found)
+    // Check field exists
+    if (fieldHeader.typeHeaderOk<fieldType>(true))
    {
-        Info<< "    - set internal values of "
+        Info<< "    Setting internal values of "
            << fieldHeader.headerClassName()
-            << ": " << fieldDesc.name()
-            << " = " << fieldValue << endl;
+            << " " << fieldName << endl;

        fieldType field(fieldHeader, mesh, false);

-        if (isNull(selectedCells) || selectedCells.size() == field.size())
+        const Type value = pTraits<Type>(fieldValueStream);
+
+        if (selectedCells.size() == field.size())
        {
-            field.primitiveFieldRef() = fieldValue;
+            field.primitiveFieldRef() = value;
        }
        else
        {
-            for (const label celli : selectedCells)
+            forAll(selectedCells, celli)
            {
-                field[celli] = fieldValue;
+                field[selectedCells[celli]] = value;
            }
        }

-        // Make boundary fields consistent - treat like zeroGradient
-        for (auto& pfld : field.boundaryFieldRef())
+        typename GeometricField<Type, fvPatchField, volMesh>::
+            Boundary& fieldBf = field.boundaryFieldRef();
+
+        forAll(field.boundaryField(), patchi)
        {
-            pfld = pfld.patchInternalField();
+            fieldBf[patchi] = fieldBf[patchi].patchInternalField();
        }

        if (!field.write())
        {
            FatalErrorInFunction
-                << "Failed writing field " << field.name() << endl;
+              << "Failed writing field " << fieldName << endl;
        }
    }
    else
    {
-        Warning
-            << "Field " << fieldDesc.name() << " not found" << endl;
+        WarningInFunction
+          << "Field " << fieldName << " not found" << endl;
+
+        // Consume value
+        (void)pTraits<Type>(fieldValueStream);
    }

    return true;
 }


-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-// Setting finite-area face fields
-template<class Type>
-bool setAreaFieldType
-(
-    const fieldDescription& fieldDesc,
-    const faMesh& mesh,
-    const labelList& selectedFaces,
-    Istream& is
-)
+class setCellField
 {
-    typedef GeometricField<Type, faPatchField, areaMesh> fieldType;

-    if (fieldDesc.type() != fieldType::typeName)
+public:
+
+    setCellField()
+    {}
+
+    autoPtr<setCellField> clone() const
    {
-        return false;
+        return autoPtr<setCellField>::New();
    }

-    // Get value from stream
-    const Type fieldValue = pTraits<Type>(is);
-
-
-    // Check the current time directory
-    IOobject fieldHeader
-    (
-        fieldDesc.name(),
-        mesh.thisDb().time().timeName(),
-        mesh.thisDb(),
-        IOobject::MUST_READ
-    );
-
-    bool found = fieldHeader.typeHeaderOk<fieldType>(true);
-
-    if (!found)
+    class iNew
    {
-        // Fallback to "constant" directory
-        fieldHeader = IOobject
-        (
-            fieldDesc.name(),
-            mesh.thisDb().time().constant(),
-            mesh.thisDb(),
-            IOobject::MUST_READ
-        );
-        found = fieldHeader.typeHeaderOk<fieldType>(true);
-    }
+        const fvMesh& mesh_;
+        const labelList& selectedCells_;

-    // Field exists
-    if (found)
-    {
-        Info<< "    - set internal values of "
-            << fieldHeader.headerClassName()
-            << ": " << fieldDesc.name()
-            << " = " << fieldValue << endl;
+    public:

-        fieldType field(fieldHeader, mesh);
+        iNew(const fvMesh& mesh, const labelList& selectedCells)
+        :
+            mesh_(mesh),
+            selectedCells_(selectedCells)
+        {}

-        if (isNull(selectedFaces) || selectedFaces.size() == field.size())
+        iNew(const fvMesh& mesh, labelList&& selectedCells)
+        :
+            mesh_(mesh),
+            selectedCells_(std::move(selectedCells))
+        {}
+
+        autoPtr<setCellField> operator()(Istream& fieldValues) const
        {
-            field.primitiveFieldRef() = fieldValue;
-        }
-        else
-        {
-            for (const label facei : selectedFaces)
+            word fieldType(fieldValues);
+
+            if
+            (
+               !(
+                    setCellFieldType<scalar>
+                        (fieldType, mesh_, selectedCells_, fieldValues)
+                 || setCellFieldType<vector>
+                        (fieldType, mesh_, selectedCells_, fieldValues)
+                 || setCellFieldType<sphericalTensor>
+                        (fieldType, mesh_, selectedCells_, fieldValues)
+                 || setCellFieldType<symmTensor>
+                        (fieldType, mesh_, selectedCells_, fieldValues)
+                 || setCellFieldType<tensor>
+                        (fieldType, mesh_, selectedCells_, fieldValues)
+                )
+            )
            {
-                field[facei] = fieldValue;
+                WarningInFunction
+                    << "field type " << fieldType << " not currently supported"
+                    << endl;
            }
+
+            return autoPtr<setCellField>::New();
        }
-
-        if (!field.write())
-        {
-            FatalErrorInFunction
-                << "Failed writing field " << field.name() << endl;
-        }
-    }
-    else
-    {
-        Warning
-            << "Field " << fieldDesc.name() << " not found" << endl;
-    }
-
-    return true;
-}
+    };
+};


-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-// Setting volume boundary fields
 template<class Type>
 bool setFaceFieldType
 (
-    const fieldDescription& fieldDesc,
+    const word& fieldTypeDesc,
    const fvMesh& mesh,
    const labelList& selectedFaces,
-    Istream& is
+    Istream& fieldValueStream
 )
 {
    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;

-    if (fieldDesc.type() != fieldType::typeName)
+    if (fieldTypeDesc != fieldType::typeName + "Value")
    {
        return false;
    }

-    // Get value from stream
-    const Type fieldValue = pTraits<Type>(is);
-
+    word fieldName(fieldValueStream);

    // Check the current time directory
    IOobject fieldHeader
    (
-        fieldDesc.name(),
-        mesh.thisDb().time().timeName(),
-        mesh.thisDb(),
+        fieldName,
+        mesh.time().timeName(),
+        mesh,
        IOobject::MUST_READ
    );

-    bool found = fieldHeader.typeHeaderOk<fieldType>(true);
-
-    if (!found)
+    // Check the "constant" directory
+    if (!fieldHeader.typeHeaderOk<fieldType>(true))
    {
-        // Fallback to "constant" directory
        fieldHeader = IOobject
        (
-            fieldDesc.name(),
-            mesh.thisDb().time().constant(),
-            mesh.thisDb(),
+            fieldName,
+            mesh.time().constant(),
+            mesh,
            IOobject::MUST_READ
        );
-        found = fieldHeader.typeHeaderOk<fieldType>(true);
    }

-    // Field exists
-    if (found)
+    // Check field exists
+    if (fieldHeader.typeHeaderOk<fieldType>(true))
    {
-        Info<< "    - set boundary values of "
+        Info<< "    Setting patchField values of "
            << fieldHeader.headerClassName()
-            << ": " << fieldDesc.name()
-            << " = " << fieldValue << endl;
+            << " " << fieldName << endl;

        fieldType field(fieldHeader, mesh);

+        const Type value = pTraits<Type>(fieldValueStream);
+
        // Create flat list of selected faces and their value.
        Field<Type> allBoundaryValues(mesh.nBoundaryFaces());
        forAll(field.boundaryField(), patchi)
@ -358,47 +260,34 @@ bool setFaceFieldType
        }

        // Override
-        unsigned hasWarned = 0;
+        bool hasWarned = false;
        labelList nChanged
        (
            returnReduce(field.boundaryField().size(), maxOp<label>()),
-            Zero
+            0
        );
-
-        for (const label facei : selectedFaces)
+        forAll(selectedFaces, i)
        {
-            const label bFacei = facei-mesh.nInternalFaces();
-
-            if (bFacei < 0)
+            label facei = selectedFaces[i];
+            if (mesh.isInternalFace(facei))
            {
-                if (!(hasWarned & 1))
+                if (!hasWarned)
                {
-                    hasWarned |= 1;
+                    hasWarned = true;
                    WarningInFunction
                        << "Ignoring internal face " << facei
                        << ". Suppressing further warnings." << endl;
                }
            }
-            else if (bFacei >= mesh.nBoundaryFaces())
-            {
-                if (!(hasWarned & 2))
-                {
-                    hasWarned |= 2;
-                    WarningInFunction
-                        << "Ignoring out-of-range face " << facei
-                        << ". Suppressing further warnings." << endl;
-                }
-            }
            else
            {
-                label patchi = mesh.boundaryMesh().patchID()[bFacei];
-
-                allBoundaryValues[bFacei] = fieldValue;
-                ++nChanged[patchi];
+                label bFacei = facei-mesh.nInternalFaces();
+                allBoundaryValues[bFacei] = value;
+                nChanged[mesh.boundaryMesh().patchID()[bFacei]]++;
            }
        }

-        Pstream::listCombineReduce(nChanged, plusEqOp<label>());
+        Pstream::listCombineAllGather(nChanged, plusEqOp<label>());

        auto& fieldBf = field.boundaryFieldRef();

@ -410,7 +299,6 @@ bool setFaceFieldType
                Info<< "    On patch "
                    << field.boundaryField()[patchi].patch().name()
                    << " set " << nChanged[patchi] << " values" << endl;
-
                fieldBf[patchi] == SubField<Type>
                (
                    allBoundaryValues,
@ -424,223 +312,86 @@ bool setFaceFieldType
        if (!field.write())
        {
            FatalErrorInFunction
-                << "Failed writing field " << field.name() << endl;
+                << "Failed writing field " << field.name() << exit(FatalError);
        }
    }
    else
    {
-        Warning
-            << "Field " << fieldDesc.name() << " not found" << endl;
+        WarningInFunction
+          << "Field " << fieldName << " not found" << endl;
+
+        // Consume value
+        (void)pTraits<Type>(fieldValueStream);
    }

    return true;
 }


-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-// Dispatcher for setting volume fields
-struct setCellField
+class setFaceField
 {
-    autoPtr<setCellField> clone() const { return nullptr; }  // placeholder

-    static bool apply
-    (
-        const fieldDescription& fieldDesc,
-        const fvMesh& m,
-        const labelList& selectedCells,
-        Istream& is
-    )
+public:
+
+    setFaceField()
+    {}
+
+    autoPtr<setFaceField> clone() const
    {
-        return
-        (
-            setCellFieldType<scalar>(fieldDesc, m, selectedCells, is)
-         || setCellFieldType<vector>(fieldDesc, m, selectedCells, is)
-         || setCellFieldType<sphericalTensor>(fieldDesc, m, selectedCells, is)
-         || setCellFieldType<symmTensor>(fieldDesc, m, selectedCells, is)
-         || setCellFieldType<tensor>(fieldDesc, m, selectedCells, is)
-        );
+        return autoPtr<setFaceField>::New();
    }

    class iNew
    {
        const fvMesh& mesh_;
-        const labelList& selected_;
-
-    public:
-
-        iNew(const fvMesh& mesh, const labelList& selectedCells)
-        :
-            mesh_(mesh),
-            selected_(selectedCells)
-        {}
-
-        autoPtr<setCellField> operator()(Istream& is) const
-        {
-            const fieldDescription fieldDesc(is);
-
-            bool ok = setCellField::apply(fieldDesc, mesh_, selected_, is);
-
-            if (!ok)
-            {
-                ok = consumeUnused(fieldDesc, is);
-
-                if (ok)
-                {
-                    // Not meant for us
-                    Info<< "Skip " << fieldDesc.type()
-                        << " for finite-volume" << nl;
-                }
-                else
-                {
-                    WarningInFunction
-                        << "Unsupported field type: "
-                        << fieldDesc.type() << endl;
-                }
-            }
-
-            return nullptr;  // Irrelevant return value
-        }
-    };
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-// Dispatcher for setting volume boundary fields
-struct setFaceField
-{
-    autoPtr<setFaceField> clone() const { return nullptr; }  // placeholder
-
-    static bool apply
-    (
-        const fieldDescription& fieldDesc,
-        const fvMesh& m,
-        const labelList& selectedFaces,
-        Istream& is
-    )
-    {
-        return
-        (
-            setFaceFieldType<scalar>(fieldDesc, m, selectedFaces, is)
-         || setFaceFieldType<vector>(fieldDesc, m, selectedFaces, is)
-         || setFaceFieldType<sphericalTensor>(fieldDesc, m, selectedFaces, is)
-         || setFaceFieldType<symmTensor>(fieldDesc, m, selectedFaces, is)
-         || setFaceFieldType<tensor>(fieldDesc, m, selectedFaces, is)
-        );
-    }
-
-    class iNew
-    {
-        const fvMesh& mesh_;
-        const labelList& selected_;
+        const labelList& selectedFaces_;

    public:

        iNew(const fvMesh& mesh, const labelList& selectedFaces)
        :
            mesh_(mesh),
-            selected_(selectedFaces)
+            selectedFaces_(selectedFaces)
        {}

-        autoPtr<setFaceField> operator()(Istream& is) const
-        {
-            const fieldDescription fieldDesc(is);
-
-            bool ok = setFaceField::apply(fieldDesc, mesh_, selected_, is);
-
-            if (!ok)
-            {
-                ok = consumeUnused(fieldDesc, is);
-
-                if (ok)
-                {
-                    // Not meant for us
-                    Info<< "Skip " << fieldDesc.type()
-                        << " for finite-volume" << nl;
-                }
-                else
-                {
-                    WarningInFunction
-                        << "Unsupported field type: "
-                        << fieldDesc.type() << endl;
-                }
-            }
-
-            return nullptr;  // Irrelevant return value
-        }
-    };
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-// Dispatcher for setting area face fields
-struct setAreaField
-{
-    autoPtr<setAreaField> clone() const { return nullptr; }  // placeholder
-
-    static bool apply
-    (
-        const fieldDescription& fieldDesc,
-        const faMesh& m,
-        const labelList& selectedFaces,
-        Istream& is
-    )
-    {
-        return
-        (
-            setAreaFieldType<scalar>(fieldDesc, m, selectedFaces, is)
-         || setAreaFieldType<vector>(fieldDesc, m, selectedFaces, is)
-         || setAreaFieldType<sphericalTensor>(fieldDesc, m, selectedFaces, is)
-         || setAreaFieldType<symmTensor>(fieldDesc, m, selectedFaces, is)
-         || setAreaFieldType<tensor>(fieldDesc, m, selectedFaces, is)
-        );
-    }
-
-    class iNew
-    {
-        const faMesh& mesh_;
-        const labelList& selected_;
-
-    public:
-
-        iNew(const faMesh& mesh, const labelList& selectedFaces)
+        iNew(const fvMesh& mesh, labelList&& selectedFaces)
        :
            mesh_(mesh),
-            selected_(selectedFaces)
+            selectedFaces_(std::move(selectedFaces))
        {}

-        autoPtr<setAreaField> operator()(Istream& is) const
+        autoPtr<setFaceField> operator()(Istream& fieldValues) const
        {
-            const fieldDescription fieldDesc(is);
+            word fieldType(fieldValues);

-            bool ok = setAreaField::apply(fieldDesc, mesh_, selected_, is);
-
-            if (!ok)
+            if
+            (
+               !(
+                    setFaceFieldType<scalar>
+                        (fieldType, mesh_, selectedFaces_, fieldValues)
+                 || setFaceFieldType<vector>
+                        (fieldType, mesh_, selectedFaces_, fieldValues)
+                 || setFaceFieldType<sphericalTensor>
+                        (fieldType, mesh_, selectedFaces_, fieldValues)
+                 || setFaceFieldType<symmTensor>
+                        (fieldType, mesh_, selectedFaces_, fieldValues)
+                 || setFaceFieldType<tensor>
+                        (fieldType, mesh_, selectedFaces_, fieldValues)
+                )
+            )
            {
-                ok = consumeUnused(fieldDesc, is);
-
-                if (ok)
-                {
-                    // Not meant for us
-                    Info<< "Skip " << fieldDesc.type()
-                        << " for finite-volume" << nl;
-                }
-                else
-                {
-                    WarningInFunction
-                        << "Unsupported field type: "
-                        << fieldDesc.type() << endl;
-                }
+                WarningInFunction
+                    << "field type " << fieldType << " not currently supported"
+                    << endl;
            }

-            return nullptr;  // Irrelevant return value
+            return autoPtr<setFaceField>::New();
        }
    };
 };


+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

 int main(int argc, char *argv[])
@ -652,28 +403,11 @@ int main(int argc, char *argv[])

    argList::addOption("dict", "file", "Alternative setFieldsDict");

-    argList::addBoolOption
-    (
-        "no-finite-area",
-        "Suppress handling of finite-area mesh/fields"
-    );
-
    #include "addRegionOption.H"
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createNamedMesh.H"

-    autoPtr<faMesh> faMeshPtr;
-
-    if (!args.found("no-finite-area"))
-    {
-        faMeshPtr = faMesh::TryNew(mesh);
-    }
-    if (faMeshPtr)
-    {
-        Info<< "Detected finite-area mesh" << nl;
-    }
-
    const word dictName("setFieldsDict");
    #include "setSystemMeshDictionaryIO.H"

@ -681,144 +415,70 @@ int main(int argc, char *argv[])

    IOdictionary setFieldsDict(dictIO);

-
-    // Note.
-    // The PtrList + iNew mechanism is used to trigger the callbacks
-    // and perform the desired actions. The contents of the PtrList
-    // itself are actually irrelevant.
-
-
-    // Default field values
+    if (setFieldsDict.found("defaultFieldValues"))
    {
-        const entry* eptr =
-            setFieldsDict.findEntry("defaultFieldValues", keyType::LITERAL);
-
-        if (eptr)
-        {
-            ITstream& is = eptr->stream();
-
-            Info<< "Setting volume field default values" << endl;
-
-            PtrList<setCellField> defaultFieldValues
-            (
-                is,
-                setCellField::iNew(mesh, labelList::null())
-            );
-
-            if (faMeshPtr)
-            {
-                const faMesh& areaMesh = faMeshPtr();
-                is.rewind();
-
-                Info<< "Setting area field default values" << endl;
-
-                PtrList<setAreaField> defaultFieldValues
-                (
-                    is,
-                    setAreaField::iNew(areaMesh, labelList::null())
-                );
-            }
-            Info<< endl;
-        }
+        Info<< "Setting field default values" << endl;
+        PtrList<setCellField> defaultFieldValues
+        (
+            setFieldsDict.lookup("defaultFieldValues"),
+            setCellField::iNew(mesh, labelList(mesh.nCells()))
+        );
+        Info<< endl;
    }


-    Info<< "Setting field region values" << nl << endl;
+    Info<< "Setting field region values" << endl;

    PtrList<entry> regions(setFieldsDict.lookup("regions"));

-    for (const entry& region : regions)
+    forAll(regions, regionI)
    {
+        const entry& region = regions[regionI];
+
        autoPtr<topoSetSource> source =
            topoSetSource::New(region.keyword(), mesh, region.dict());

        if (source().setType() == topoSetSource::CELLSET_SOURCE)
        {
-            cellSet subset
+            cellSet selectedCellSet
            (
                mesh,
                "cellSet",
                mesh.nCells()/10+1  // Reasonable size estimate.
            );

-            source->applyToSet(topoSetSource::NEW, subset);
-
-            labelList selectedCells(subset.sortedToc());
-
-            Info<< "    Selected "
-                << returnReduce(selectedCells.size(), sumOp<label>())
-                << '/'
-                << returnReduce(mesh.nCells(), sumOp<label>())
-                << " cells" << nl;
-
-            ITstream& is = region.dict().lookup("fieldValues");
+            source->applyToSet
+            (
+                topoSetSource::NEW,
+                selectedCellSet
+            );

            PtrList<setCellField> fieldValues
            (
-                is,
-                setCellField::iNew(mesh, selectedCells)
+                region.dict().lookup("fieldValues"),
+                setCellField::iNew(mesh, selectedCellSet.sortedToc())
            );
        }
        else if (source().setType() == topoSetSource::FACESET_SOURCE)
        {
-            faceSet subset
+            faceSet selectedFaceSet
            (
                mesh,
                "faceSet",
                mesh.nBoundaryFaces()/10+1
            );

-            source->applyToSet(topoSetSource::NEW, subset);
-
-            labelList selectedFaces(subset.sortedToc());
-
-            Info<< "    Selected " << selectedFaces.size()
-                << " faces" << nl;
-
-            ITstream& is = region.dict().lookup("fieldValues");
+            source->applyToSet
+            (
+                topoSetSource::NEW,
+                selectedFaceSet
+            );

            PtrList<setFaceField> fieldValues
            (
-                is,
-                setFaceField::iNew(mesh, selectedFaces)
+                region.dict().lookup("fieldValues"),
+                setFaceField::iNew(mesh, selectedFaceSet.sortedToc())
            );
-
-            if (faMeshPtr)
-            {
-                const faMesh& areaMesh = faMeshPtr();
-
-                const labelUList& faceLabels = areaMesh.faceLabels();
-
-                // Transcribe from mesh faces to finite-area addressing
-                labelList areaFaces(faceLabels.size());
-
-                label nUsed = 0;
-                forAll(faceLabels, facei)
-                {
-                    const label meshFacei = faceLabels[facei];
-
-                    if (subset.test(meshFacei))
-                    {
-                        areaFaces[nUsed] = facei;
-                        ++nUsed;
-                    }
-                }
-                areaFaces.resize(nUsed);
-
-                Info<< "    Selected "
-                    << returnReduce(areaFaces.size(), sumOp<label>())
-                    << '/'
-                    << returnReduce(faceLabels.size(), sumOp<label>())
-                    << " area faces" << nl;
-
-                is.rewind();
-
-                PtrList<setAreaField> fieldValues
-                (
-                    is,
-                    setAreaField::iNew(areaMesh, areaFaces)
-                );
-            }
        }
    }

--- a/applications/utilities/preProcessing/setTurbulenceFields/setTurbulenceFields.C
+++ b/applications/utilities/preProcessing/setTurbulenceFields/setTurbulenceFields.C
@ -109,7 +109,6 @@ Note
 #include "singlePhaseTransportModel.H"
 #include "turbulentTransportModel.H"
 #include "turbulentFluidThermoModel.H"
-#include "processorFvPatchField.H"
 #include "wallFvPatch.H"
 #include "fixedValueFvPatchFields.H"

@ -123,40 +122,6 @@ void InfoField(const word& fldName)
 }


-template<class Type>
-void correctProcessorPatches
-(
-    GeometricField<Type, fvPatchField, volMesh>& vf
-)
-{
-    if (!Pstream::parRun())
-    {
-        return;
-    }
-
-    // Not possible to use correctBoundaryConditions on fields as they may
-    // use local info as opposed to the constraint values employed here,
-    // but still need to update processor patches
-    auto& bf = vf.boundaryFieldRef();
-
-    forAll(bf, patchi)
-    {
-        if (isA<processorFvPatchField<Type>>(bf[patchi]))
-        {
-            bf[patchi].initEvaluate();
-        }
-    }
-
-    forAll(bf, patchi)
-    {
-        if (isA<processorFvPatchField<Type>>(bf[patchi]))
-        {
-            bf[patchi].evaluate();
-        }
-    }
-}
-
-
 IOobject createIOobject
 (
    const fvMesh& mesh,
@ -482,26 +447,22 @@ int main(int argc, char *argv[])
        // (M:Eq. 9)
        const dimensionedScalar maxU(dimVelocity, SMALL);
        U *= min(scalar(1), fRei*uTau/max(mag(U), maxU));
-        correctProcessorPatches<vector>(U);
    }

    if (tepsilon.valid())
    {
        tepsilon.ref() = epsilon;
-        correctProcessorPatches<scalar>(tepsilon.ref());
    }

    if (tk.valid())
    {
        tk.ref() = k;
-        correctProcessorPatches<scalar>(tk.ref());
    }

    if (tomega.valid())
    {
        const dimensionedScalar k0(sqr(dimLength/dimTime), SMALL);
        tomega.ref() = Cmu*epsilon/(k + k0);
-        correctProcessorPatches<scalar>(tomega.ref());
    }

    if (tR.valid())
@ -514,7 +475,6 @@ int main(int argc, char *argv[])
        {
            R[celli] = Rdiag[celli];
        }
-        correctProcessorPatches<symmTensor>(R);
    }


--- a/applications/utilities/preProcessing/viewFactorsGen/Allwmake
+++ b/applications/utilities/preProcessing/viewFactorsGen/Allwmake
@ -1,19 +0,0 @@
-#!/bin/sh
-cd "${0%/*}" || exit                                # Run from this directory
-. ${WM_PROJECT_DIR:?}/wmake/scripts/AllwmakeParseArguments # (error catching)
-. ${WM_PROJECT_DIR:?}/wmake/scripts/have_cgal
-
-#------------------------------------------------------------------------------
-unset COMP_FLAGS LINK_FLAGS
-
-if have_cgal
-then
-    echo "    found CGAL -- enabling CGAL support."
-else
-    echo "    did not find CGAL -- disabling CGAL functionality"
-    export COMP_FLAGS="-DNO_CGAL"
-fi
-
-wmake $targetType
-
-#------------------------------------------------------------------------------
--- a/applications/utilities/preProcessing/viewFactorsGen/Make/options
+++ b/applications/utilities/preProcessing/viewFactorsGen/Make/options
@ -1,18 +1,10 @@
-include $(GENERAL_RULES)/cgal-header-only
-
 EXE_INC = \
-    -Wno-old-style-cast \
-    $(COMP_FLAGS) \
-    ${CGAL_INC} \
-    -DCGAL_HEADER_ONLY \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/surfMesh/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/parallel/distributed/lnInclude
-
+    -I$(LIB_SRC)/parallel/distributed/lnInclude \

 EXE_LIBS = \
-    /* ${CGAL_LIBS} */ \
    -lfiniteVolume \
    -lsurfMesh \
    -lmeshTools \
--- a/applications/utilities/preProcessing/viewFactorsGen/searchingEngine_CGAL.H
+++ b/applications/utilities/preProcessing/viewFactorsGen/searchingEngine_CGAL.H
@ -1,84 +0,0 @@
-Random rndGen(653213);
-
-// Determine mesh bounding boxes:
-List<treeBoundBox> meshBb
-(
-    1,
-    treeBoundBox
-    (
-        boundBox(coarseMesh.points(), false)
-    ).extend(rndGen, 1e-3)
-);
-
-// Dummy bounds dictionary
-dictionary dict;
-dict.add("bounds", meshBb);
-dict.add
-(
-    "distributionType",
-    distributedTriSurfaceMesh::distributionTypeNames_
-    [
-        distributedTriSurfaceMesh::FROZEN
-    ]
-);
-dict.add("mergeDistance", SMALL);
-
-labelList triSurfaceToAgglom(5*nFineFaces);
-
-triSurface localSurface = triangulate
-(
-    patches,
-    includePatches,
-    finalAgglom,
-    triSurfaceToAgglom,
-    globalNumbering,
-    coarsePatches
-);
-
-
-// CGAL surface
-
-distributedTriSurfaceMesh surfacesMesh
-(
-    IOobject
-    (
-        "wallSurface.stl",
-        runTime.constant(),     // directory
-        "triSurface",           // instance
-        runTime,                // registry
-        IOobject::NO_READ,
-        IOobject::NO_WRITE
-    ),
-    localSurface,
-    dict
-);
-
-triSurfaceToAgglom.resize(surfacesMesh.size());
-
-surfacesMesh.setField(triSurfaceToAgglom);
-
-const pointField& pts = surfacesMesh.localPoints();
-
-std::list<Triangle> triangles;
-
-for (const auto& triLocal : surfacesMesh.localFaces())
-{
-    point p1l = pts[triLocal[0]];
-    point p2l = pts[triLocal[1]];
-    point p3l = pts[triLocal[2]];
-
-    Point p1(p1l[0], p1l[1], p1l[2]);
-    Point p2(p2l[0], p2l[1], p2l[2]);
-    Point p3(p3l[0], p3l[1], p3l[2]);
-
-    Triangle tri(p1, p2, p3);
-
-    if (tri.is_degenerate())
-    {
-        std::cout << tri << std::endl;
-    }
-    triangles.push_back(tri);
-}
-
-// constructs AABB tree
-Tree tree(triangles.begin(), triangles.end());
--- a/applications/utilities/preProcessing/viewFactorsGen/shootRays.H
+++ b/applications/utilities/preProcessing/viewFactorsGen/shootRays.H
@ -1,15 +1,6 @@
 // All rays expressed as start face (local) index end end face (global)
 // Pre-size by assuming a certain percentage is visible.

-
-if (Pstream::master())
-{
-    Info<< "\nShooting rays using distributedTriSurfaceMesh (no CGAL support)"
-        << endl;
-}
-
-
-
 // Maximum length for dynamicList
 const label maxDynListLength
 (
--- a/applications/utilities/preProcessing/viewFactorsGen/shootRays_CGAL.H
+++ b/applications/utilities/preProcessing/viewFactorsGen/shootRays_CGAL.H
@ -1,131 +0,0 @@
-// Maximum length for dynamicList
-const label maxDynListLength
-(
-    viewFactorDict.getOrDefault<label>("maxDynListLength", 1000000000)
-);
-
-for (const int proci : Pstream::allProcs())
-{
-    std::vector<Point> start;
-    start.reserve(coarseMesh.nFaces());
-
-    std::vector<Point> end(start.size());
-    end.reserve(start.size());
-
-    DynamicList<label> startIndex(start.size());
-    DynamicList<label> endIndex(start.size());
-
-    DynamicList<label> startAgg(start.size());
-    DynamicList<label> endAgg(start.size());
-
-    const pointField& myFc = remoteCoarseCf[Pstream::myProcNo()];
-    const vectorField& myArea = remoteCoarseSf[Pstream::myProcNo()];
-    const labelField& myAgg = remoteCoarseAgg[Pstream::myProcNo()];
-
-    const pointField& remoteArea = remoteCoarseSf[proci];
-    const pointField& remoteFc = remoteCoarseCf[proci];
-    const labelField& remoteAgg = remoteCoarseAgg[proci];
-
-    label i = 0;
-    label j = 0;
-    do
-    {
-        for (; i < myFc.size(); i++)
-        {
-            const point& fc = myFc[i];
-            const vector& fA = myArea[i];
-            const label& fAgg = myAgg[i];
-
-            for (; j < remoteFc.size(); j++)
-            {
-                if (proci != Pstream::myProcNo() || i != j)
-                {
-                    const point& remFc = remoteFc[j];
-                    const vector& remA = remoteArea[j];
-                    const label& remAgg = remoteAgg[j];
-
-                    const vector d(remFc - fc);
-
-
-                    const vector nd = d/mag(d);
-                    const vector nfA = fA/mag(fA);
-                    const vector nremA = remA/mag(remA);
-
-                    if (((nd & nfA) < 0) && ((nd & nremA) > 0))
-                    {
-                        vector direction(d[0], d[1], d[2]);
-                        vector s(fc[0], fc[1], fc[2]);
-                        vector rayEnd(s + (1-intTol)*direction);
-                        end.push_back(Point(rayEnd[0], rayEnd[1], rayEnd[2]));
-
-                        s += vector(intTol*d[0], intTol*d[1], intTol*d[2]);
-
-                        start.push_back(Point(s[0], s[1], s[2]));
-                        startIndex.append(i);
-                        if (useAgglomeration)
-                        {
-                            startAgg.append
-                            (
-                                globalNumbering.toGlobal(proci, fAgg)
-                            );
-                            endAgg.append
-                            (
-                                globalNumbering.toGlobal(proci, remAgg)
-                            );
-                        }
-
-                        label globalI = globalNumbering.toGlobal(proci, j);
-                        endIndex.append(globalI);
-
-                        if (startIndex.size() > maxDynListLength)
-                        {
-                            FatalErrorInFunction
-                                << "Dynamic list need from capacity."
-                                << "Actual size maxDynListLength : "
-                                <<  maxDynListLength
-                                << abort(FatalError);
-                        }
-                    }
-                }
-            }
-
-            if (j == remoteFc.size())
-            {
-                j = 0;
-            }
-        }
-
-    } while (i < myFc.size());
-
-    label totalRays(startIndex.size());
-    reduce(totalRays, sumOp<scalar>());
-
-    if (debug)
-    {
-        Pout<< "Number of rays :" << totalRays << endl;
-    }
-
-    for (unsigned long rayI = 0; rayI < start.size(); ++rayI)
-    {
-        Segment ray(start[rayI], end[rayI]);
-        Segment_intersection intersects = tree.any_intersection(ray);
-
-        if (!intersects)
-        {
-            rayStartFace.append(startIndex[rayI]);
-            rayEndFace.append(endIndex[rayI]);
-        }
-    }
-    start.clear();
-
-    if (debug)
-    {
-        Pout << "hits : " <<  rayStartFace.size()<< endl;
-    }
-
-    startIndex.clear();
-    end.clear();
-    endIndex.clear();
-    startAgg.clear();
-    endAgg.clear();
-}
--- a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
+++ b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
@ -31,34 +31,14 @@ Group
    grpPreProcessingUtilities

 Description
-    This view factors generation application uses a combined approach of
-    double area integral (2AI) and double linear integral (2LI). 2AI is used
-    when the two surfaces are 'far' apart and 2LI when they are 'close'.
-    2LI is integrated along edges using Gaussian quadrature.
-    The distance between faces is calculating a ratio between averaged areas
-    and the distance between face centres.
+    View factors are calculated based on a face agglomeration array
+    (finalAgglom generated by faceAgglomerate utility).

-    The input from viewFactorsDict are:
+    Each view factor between the agglomerated faces i and j (Fij) is calculated
+    using a double integral of the sub-areas composing the agglomeration.

-        GaussQuadTol              0.1;      // GaussQuad  error
-        distTol                   8;        // R/Average(rm)
-        alpha                     0.22;     // Use for common edges for 2LI
-
-
-    For debugging purposes, the following entries can be set in viewFactorsDict:
-
-        writeViewFactorMatrix     true;
-        writeFacesAgglomeration   false;
-        dumpRays		          false;
-
-
-        writeViewFactorMatrix   writes the sum of the VF on each face.
-        writeFacesAgglomeration writes the agglomeration
-        dumpRays                dumps rays
-
-
-    The participating patches in the VF calculation have to be in the
-    'viewFactorWall' patch group (in the polyMesh/boundary file), e.g.
+    The patches involved in the view factor calculation are taken from the
+    boundary file and should be part on the group viewFactorWall. ie.:

    floor
    {
@ -68,9 +48,6 @@ Description
        startFace       3100;
    }

-    Compile with -DNO_CGAL only if no CGAL present - CGAL AABB tree performs
-    better than the built-in octree.
-
 \*---------------------------------------------------------------------------*/

 #include "argList.H"
@ -79,52 +56,21 @@ Description
 #include "surfaceFields.H"
 #include "distributedTriSurfaceMesh.H"
 #include "meshTools.H"
-#include "constants.H"

 #include "indirectPrimitivePatch.H"
 #include "DynamicField.H"
+#include "unitConversion.H"

 #include "scalarMatrices.H"
 #include "labelListIOList.H"
 #include "scalarListIOList.H"

 #include "singleCellFvMesh.H"
+
 #include "IOmapDistribute.H"

-#ifndef NO_CGAL
-
-// Silence boost bind deprecation warnings (before CGAL-5.2.1)
-#include "CGAL/version.h"
-#if defined(CGAL_VERSION_NR) && (CGAL_VERSION_NR < 1050211000)
-#define BOOST_BIND_GLOBAL_PLACEHOLDERS
-#endif
-
-#include <CGAL/Simple_cartesian.h>
-#include <CGAL/AABB_tree.h>
-#include <CGAL/AABB_traits.h>
-#include <CGAL/AABB_triangle_primitive.h>
-#include <CGAL/Surface_mesh.h>
-
-typedef CGAL::Simple_cartesian<double> K;
-typedef K::Point_3 Point;
-typedef K::Direction_3 Vector3C;
-typedef K::Triangle_3 Triangle;
-typedef K::Segment_3 Segment;
-
-typedef std::list<Triangle>::iterator Iterator;
-typedef CGAL::AABB_triangle_primitive<K, Iterator> Primitive;
-typedef CGAL::AABB_traits<K, Primitive> AABB_triangle_traits;
-typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
-typedef boost::optional
-<
-    Tree::Intersection_and_primitive_id<Segment>::Type
-> Segment_intersection;
-
-#endif // NO_CGAL
-
 using namespace Foam;
-using namespace Foam::constant;
-using namespace Foam::constant::mathematical;
+

 triSurface triangulate
 (
@ -181,56 +127,20 @@ triSurface triangulate
        newPatchI++;
    }

-    //triSurfaceToAgglom.resize(localTriFaceI-1);
+    //striSurfaceToAgglom.resize(localTriFaceI-1);

    triangles.shrink();
-    triSurface surface(triangles, mesh.points());
-    surface.compactPoints();

+    // Create globally numbered tri surface
+    triSurface rawSurface(triangles, mesh.points());

-#ifndef NO_CGAL
-
-    // CGAL : every processor has whole surface
-
-    globalIndex globalFaceIdx(surface.size(), globalIndex::gatherOnly());
-    globalIndex globalPointIdx
+    // Create locally numbered tri surface
+    triSurface surface
    (
-        surface.points().size(), globalIndex::gatherOnly()
+        rawSurface.localFaces(),
+        rawSurface.localPoints()
    );

-    List<labelledTri> globalSurfaceTris(globalFaceIdx.gather(surface));
-    pointField globalSurfacePoints(globalPointIdx.gather(surface.points()));
-
-    //label offset = 0;
-    for (const label proci : globalPointIdx.allProcs())
-    {
-        const label offset = globalPointIdx.localStart(proci);
-
-        if (offset)
-        {
-            for
-            (
-                labelledTri& tri
-             :  globalSurfaceTris.slice(globalFaceIdx.range(proci))
-            )
-            {
-                tri[0] += offset;
-                tri[1] += offset;
-                tri[2] += offset;
-            }
-        }
-    }
-
-    surface =
-        triSurface
-        (
-            std::move(globalSurfaceTris),
-            std::move(globalSurfacePoints)
-        );
-
-    Pstream::broadcast(surface);
-#endif
-
    // Add patch names to surface
    surface.patches().setSize(newPatchI);

@ -283,7 +193,7 @@ void writeRays
 }


-scalar calculateViewFactorFij2AI
+scalar calculateViewFactorFij
 (
    const vector& i,
    const vector& j,
@ -340,74 +250,6 @@ void insertMatrixElements
 }


-scalar GaussQuad
-(
-
-    const scalarList& w,
-    const scalarList& p,
-    const scalar& magSi,
-    const scalar& magSj,
-    const vector& di,
-    const vector& dj,
-    const vector& ci,
-    const vector& cj,
-    const scalar cosij,
-    const scalar alpha,
-    label gi
-)
-{
-    scalar dIntFij = 0;
-    if (gi == 0)
-    {
-        vector r(ci - cj);
-        if (mag(r) < SMALL)
-        {
-            r = (alpha*magSi)*di;
-        }
-        dIntFij = max(cosij*Foam::log(r&r)*magSi*magSj, 0);
-    }
-    else
-    {
-        List<vector> pi(w.size());
-        forAll (pi, i)
-        {
-            pi[i] = ci + p[i]*(magSi/2)*di;
-        }
-
-        List<vector> pj(w.size());
-        forAll (pj, i)
-        {
-            pj[i] = cj + p[i]*(magSj/2)*dj;
-        }
-
-        forAll (w, i)
-        {
-            forAll (w, j)
-            {
-                vector r(pi[i] - pj[j]);
-
-                if (mag(r) < SMALL)
-                {
-                    r = (alpha*magSi)*di;
-                     dIntFij +=
-                        cosij*w[i]*w[j]*Foam::log(r&r);
-                }
-                else
-                {
-                    dIntFij +=
-                        cosij*w[i]*w[j]*Foam::log(r&r);
-                }
-
-            }
-        }
-
-        dIntFij *= (magSi/2) * (magSj/2);
-
-    }
-    return dIntFij;
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

 int main(int argc, char *argv[])
@ -431,7 +273,7 @@ int main(int argc, char *argv[])
            "viewFactorsDict",
            runTime.constant(),
            mesh,
-            IOobject::MUST_READ,
+            IOobject::MUST_READ_IF_MODIFIED,
            IOobject::NO_WRITE
       )
    );
@ -446,23 +288,6 @@ int main(int argc, char *argv[])

    const label debug = viewFactorDict.getOrDefault<label>("debug", 0);

-    const scalar GaussQuadTol =
-        viewFactorDict.getOrDefault<scalar>("GaussQuadTol", 0.01);
-
-    const scalar distTol =
-         viewFactorDict.getOrDefault<scalar>("distTol", 8);
-
-    const scalar alpha =
-         viewFactorDict.getOrDefault<scalar>("alpha", 0.21);
-
-    const scalar intTol =
-         viewFactorDict.getOrDefault<scalar>("intTol", 1e-2);
-
-    bool useAgglomeration(true);
-
-    const polyBoundaryMesh& patches = mesh.boundaryMesh();
-    const labelList viewFactorsPatches(patches.indices(viewFactorWall));
-
    // Read agglomeration map
    labelListIOList finalAgglom
    (
@ -471,22 +296,12 @@ int main(int argc, char *argv[])
            "finalAgglom",
            mesh.facesInstance(),
            mesh,
-            IOobject::READ_IF_PRESENT,
+            IOobject::MUST_READ,
            IOobject::NO_WRITE,
            false
        )
    );

-    if (!finalAgglom.typeHeaderOk<labelListIOList>())
-    {
-        finalAgglom.setSize(patches.size());
-        for (label patchi=0;  patchi < patches.size(); patchi++)
-        {
-            finalAgglom[patchi] = identity(patches[patchi].size());
-        }
-        useAgglomeration = false;
-    }
-
    // Create the coarse mesh  using agglomeration
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

@ -521,8 +336,10 @@ int main(int argc, char *argv[])
    label nCoarseFaces = 0;      //total number of coarse faces
    label nFineFaces = 0;        //total number of fine faces

+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
    const polyBoundaryMesh& coarsePatches = coarseMesh.boundaryMesh();

+    labelList viewFactorsPatches(patches.indices(viewFactorWall));
    for (const label patchi : viewFactorsPatches)
    {
        nCoarseFaces += coarsePatches[patchi].size();
@ -648,13 +465,7 @@ int main(int argc, char *argv[])

    // Set up searching engine for obstacles
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-#ifdef NO_CGAL
-    // Using octree
    #include "searchingEngine.H"
-#else
-    // Using CGAL aabbtree (faster, more robust)
-    #include "searchingEngine_CGAL.H"
-#endif

    // Determine rays between coarse face centres
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -662,15 +473,11 @@ int main(int argc, char *argv[])

    DynamicList<label> rayEndFace(rayStartFace.size());

+
    // Return rayStartFace in local index and rayEndFace in global index
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-#ifdef NO_CGAL
-    // Using octree, distributedTriSurfaceMesh
+
    #include "shootRays.H"
-#else
-    // Using CGAL aabbtree (faster, more robust)
-    #include "shootRays_CGAL.H"
-#endif

    // Calculate number of visible faces from local index
    labelList nVisibleFaceFaces(nCoarseFaces, Zero);
@ -711,9 +518,9 @@ int main(int argc, char *argv[])
        visibleFaceFaces[faceI][nVisibleFaceFaces[faceI]++] = compactI;
    }

+
    // Construct data in compact addressing
-    // (2AA) need coarse (Ai), fine Sf (dAi) and fine Cf(r) to calculate Fij
-    // (2LI) need edges (li)
+    // I need coarse Sf (Ai), fine Sf (dAi) and fine Cf(r) to calculate Fij
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    pointField compactCoarseCf(map.constructSize(), Zero);
@ -721,16 +528,12 @@ int main(int argc, char *argv[])
    List<List<point>> compactFineSf(map.constructSize());
    List<List<point>> compactFineCf(map.constructSize());

-    DynamicList<List<point>> compactPoints(map.constructSize());
-
    DynamicList<label> compactPatchId(map.constructSize());

    // Insert my coarse local values
    SubList<point>(compactCoarseSf, nCoarseFaces) = localCoarseSf;
    SubList<point>(compactCoarseCf, nCoarseFaces) = localCoarseCf;

-    const faceList& faces = mesh.faces();
-
    // Insert my fine local values
    label compactI = 0;
    forAll(viewFactorsPatches, i)
@ -745,21 +548,11 @@ int main(int argc, char *argv[])
            const labelList& coarsePatchFace =
                coarseMesh.patchFaceMap()[patchID];

-            const polyPatch& pp = patches[patchID];
-
            forAll(coarseToFine, coarseI)
            {
                compactPatchId.append(patchID);
                List<point>& fineCf = compactFineCf[compactI];
-                List<point>& fineSf = compactFineSf[compactI];
-
-                label startFace = pp.start();
-
-                const vectorField locPoints
-                (
-                    mesh.points(),
-                    faces[coarseI + startFace]
-                );
+                List<point>& fineSf = compactFineSf[compactI++];

                const label coarseFaceI = coarsePatchFace[coarseI];
                const labelList& fineFaces = coarseToFine[coarseFaceI];
@ -767,8 +560,6 @@ int main(int argc, char *argv[])
                fineCf.setSize(fineFaces.size());
                fineSf.setSize(fineFaces.size());

-                compactPoints.append(locPoints);
-
                fineCf = UIndirectList<point>
                (
                    mesh.Cf().boundaryField()[patchID],
@ -779,25 +570,19 @@ int main(int argc, char *argv[])
                    mesh.Sf().boundaryField()[patchID],
                    coarseToFine[coarseFaceI]
                );
-
-                compactI++;
            }
        }
    }

-    if (Pstream::master() && debug)
-    {
-        Info<< "map distribute..."  << endl;
-    }
-
    // Do all swapping
    map.distribute(compactCoarseSf);
    map.distribute(compactCoarseCf);
    map.distribute(compactFineCf);
    map.distribute(compactFineSf);
-    map.distribute(compactPoints);
+
    map.distribute(compactPatchId);

+
    // Plot all rays between visible faces.
    if (dumpRays)
    {
@ -813,7 +598,8 @@ int main(int argc, char *argv[])

    // Fill local view factor matrix
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    scalarListIOList F2LI
+
+    scalarListIOList F
    (
        IOobject
        (
@ -844,61 +630,6 @@ int main(int argc, char *argv[])
        Info<< "\nCalculating view factors..." << endl;
    }

-    FixedList<scalarList, 5> GaussPoints;
-    GaussPoints[0].setSize(1);
-    GaussPoints[0] = 0;
-
-    GaussPoints[1].setSize(2);
-    GaussPoints[1][0] =  1/std::sqrt(3);
-    GaussPoints[1][1] = -1/std::sqrt(3);
-
-    GaussPoints[2].setSize(3);
-    GaussPoints[2][0] =  0;
-    GaussPoints[2][1] =  std::sqrt(3.0/5.0);
-    GaussPoints[2][2] = -std::sqrt(3.0/5.0);
-
-    GaussPoints[3].setSize(4);
-    GaussPoints[3][0] = std::sqrt(3.0/7.0 - (2.0/7.0)*std::sqrt(6.0/5.0));
-    GaussPoints[3][1] = -GaussPoints[3][0];
-    GaussPoints[3][2] = std::sqrt(3.0/7.0 + (2.0/7.0)*std::sqrt(6.0/5.0));
-    GaussPoints[3][3] = -GaussPoints[3][2];
-
-    GaussPoints[4].setSize(5);
-    GaussPoints[4][0] =  0;
-    GaussPoints[4][1] = (1.0/3.0)*std::sqrt(5.0 - 2.0*std::sqrt(10.0/7.0));
-    GaussPoints[4][2] = -GaussPoints[4][1];
-    GaussPoints[4][3] = (1.0/3.0)*std::sqrt(5.0 + 2.0*std::sqrt(10.0/7.0));
-    GaussPoints[4][4] = -GaussPoints[4][3];
-
-
-    List<scalarList> GaussWeights(5);
-    GaussWeights[0].setSize(1);
-    GaussWeights[0] = 2;
-
-    GaussWeights[1].setSize(2);
-    GaussWeights[1][0] =  1;
-    GaussWeights[1][1] =  1;
-
-    GaussWeights[2].setSize(3);
-    GaussWeights[2][0] =  8.0/9.0;
-    GaussWeights[2][1] =  5.0/9.0;
-    GaussWeights[2][2] =  5.0/9.0;
-
-    GaussWeights[3].setSize(4);
-    GaussWeights[3][0] =  (18.0 + std::sqrt(30))/36.0;
-    GaussWeights[3][1] =  (18.0 + std::sqrt(30))/36.0;
-    GaussWeights[3][2] =  (18.0 - std::sqrt(30))/36.0;
-    GaussWeights[3][3] =  (18.0 - std::sqrt(30))/36.0;
-
-    GaussWeights[4].setSize(5);
-    GaussWeights[4][0] =  128.0/225.0;
-    GaussWeights[4][1] =  (322.0 + 13.0*std::sqrt(70))/900.0;
-    GaussWeights[4][2] =  (322.0 + 13.0*std::sqrt(70))/900.0;
-    GaussWeights[4][3] =  (322.0 - 13.0*std::sqrt(70))/900.0;
-    GaussWeights[4][4] =  (322.0 - 13.0*std::sqrt(70))/900.0;
-
-    const label maxQuadOrder = 5;
-
    if (mesh.nSolutionD() == 3)
    {
        forAll(localCoarseSf, coarseFaceI)
@ -907,212 +638,119 @@ int main(int argc, char *argv[])
            const vector Ai = sum(localFineSf);
            const List<point>& localFineCf = compactFineCf[coarseFaceI];
            const label fromPatchId = compactPatchId[coarseFaceI];
-
-            const List<point>& lPoints = compactPoints[coarseFaceI];
-
            patchArea[fromPatchId] += mag(Ai);

            const labelList& visCoarseFaces = visibleFaceFaces[coarseFaceI];

            forAll(visCoarseFaces, visCoarseFaceI)
            {
-                //F2AI[coarseFaceI].setSize(visCoarseFaces.size());
-                F2LI[coarseFaceI].setSize(visCoarseFaces.size());
+                F[coarseFaceI].setSize(visCoarseFaces.size());
                label compactJ = visCoarseFaces[visCoarseFaceI];
                const List<point>& remoteFineSj = compactFineSf[compactJ];
                const List<point>& remoteFineCj = compactFineCf[compactJ];

-                const List<point>& rPointsCj = compactPoints[compactJ];
-
                const label toPatchId = compactPatchId[compactJ];

-                bool far(false);
-                // Relative distance
+                scalar Fij = 0;
                forAll(localFineSf, i)
                {
-                    const scalar dAi =
-                        Foam::sqrt
-                        (
-                            mag(localFineSf[i])/constant::mathematical::pi
-                        );
+                    const vector& dAi = localFineSf[i];
                    const vector& dCi = localFineCf[i];

                    forAll(remoteFineSj, j)
                    {
-                        const scalar dAj =
-                            Foam::sqrt
-                            (
-                                mag(remoteFineSj[j])/constant::mathematical::pi
-                            );
+                        const vector& dAj = remoteFineSj[j];
                        const vector& dCj = remoteFineCj[j];

-                        const scalar dist = mag(dCi - dCj)/((dAi + dAj)/2);
-
-                        if (dist > distTol)
-                        {
-                            far = true;
-                        }
-                    }
-                }
-
-                if (far)
-                {
-                    // 2AI method
-                    scalar F2AIij = 0;
-
-                    forAll(localFineSf, i)
-                    {
-                        const vector& dAi = localFineSf[i];
-                        const vector& dCi = localFineCf[i];
-
-                        forAll(remoteFineSj, j)
-                        {
-                            const vector& dAj = remoteFineSj[j];
-                            const vector& dCj = remoteFineCj[j];
-
-                            scalar dIntFij = calculateViewFactorFij2AI
-                            (
-                                dCi,
-                                dCj,
-                                dAi,
-                                dAj
-                            );
-
-                            F2AIij += dIntFij;
-                        }
-                    }
-                    F2LI[coarseFaceI][visCoarseFaceI] = F2AIij/mag(Ai);
-                }
-                else
-                {
-                    // 2LI method
-                    label nLocal = lPoints.size();
-                    label nRemote = rPointsCj.size();
-
-                    // Using sub-divisions (quadrature)
-                    scalar oldEToeInt = 0;
-                    for (label gi=0; gi < maxQuadOrder; gi++)
-                    {
-                        scalar F2LIij = 0;
-                        for(label i=0; i<nLocal; i++)
-                        {
-                            vector si;
-                            vector ci;
-
-                            vector sj;
-                            vector cj;
-
-                            if (i == 0)
-                            {
-                                si = lPoints[i] - lPoints[nLocal-1];
-                                ci = (lPoints[i] + lPoints[nLocal-1])/2;
-                            }
-                            else
-                            {
-                                si = lPoints[i] - lPoints[i-1];
-                                ci = (lPoints[i] + lPoints[i-1])/2;
-                            }
-
-                            for(label j=0; j<nRemote; j++)
-                            {
-                                if (j == 0)
-                                {
-                                    sj = rPointsCj[j]-rPointsCj[nRemote-1];
-                                    cj = (rPointsCj[j]+rPointsCj[nRemote-1])/2;
-                                }
-                                else
-                                {
-                                    sj = rPointsCj[j] - rPointsCj[j-1];
-                                    cj = (rPointsCj[j] + rPointsCj[j-1])/2;
-                                }
-
-
-                                scalar magSi = mag(si);
-                                scalar magSj = mag(sj);
-                                scalar cosij = (si & sj)/(magSi * magSj);
-
-                                vector di = si/magSi;
-                                vector dj = sj/magSj;
-
-                                label quadOrder = gi;
-                                const vector r(ci - cj);
-                                // Common edges use n = 0
-                                if (mag(r) < SMALL)
-                                {
-                                    quadOrder = 0;
-                                }
-
-                                scalar dIntFij =
-                                    GaussQuad
-                                    (
-                                        GaussWeights[gi],
-                                        GaussPoints[gi],
-                                        magSi,
-                                        magSj,
-                                        di,
-                                        dj,
-                                        ci,
-                                        cj,
-                                        cosij,
-                                        alpha,
-                                        quadOrder
-                                    );
-
-                                F2LIij += dIntFij;
-                            }
-                        }
-
-                        const scalar err
+                        scalar dIntFij = calculateViewFactorFij
                        (
-                            mag(F2LIij) > ROOTVSMALL
-                          ? (F2LIij-oldEToeInt)/F2LIij
-                          : 0
+                            dCi,
+                            dCj,
+                            dAi,
+                            dAj
                        );

-                        if
-                        (
-                            (mag(err) < GaussQuadTol && gi > 0)
-                         || gi == maxQuadOrder-1
-                        )
-                        {
-                            F2LI[coarseFaceI][visCoarseFaceI] =
-                                F2LIij/mag(Ai)/4/constant::mathematical::pi;
-                            break;
-                        }
-                        else
-                        {
-                            oldEToeInt = F2LIij;
-                        }
+                        Fij += dIntFij;
                    }
                }
-
-                sumViewFactorPatch[fromPatchId][toPatchId] +=
-                    F2LI[coarseFaceI][visCoarseFaceI]*mag(Ai);
+                F[coarseFaceI][visCoarseFaceI] = Fij/mag(Ai);
+                sumViewFactorPatch[fromPatchId][toPatchId] += Fij;
            }
        }
    }
-    else
+    else if (mesh.nSolutionD() == 2)
    {
-         FatalErrorInFunction
-            << " View factors are not available in 2D "
-            << exit(FatalError);
+        const boundBox& box = mesh.bounds();
+        const Vector<label>& dirs = mesh.geometricD();
+        vector emptyDir = Zero;
+        forAll(dirs, i)
+        {
+            if (dirs[i] == -1)
+            {
+                emptyDir[i] = 1.0;
+            }
+        }
+
+        scalar wideBy2 = (box.span() & emptyDir)*2.0;
+
+        forAll(localCoarseSf, coarseFaceI)
+        {
+            const vector& Ai = localCoarseSf[coarseFaceI];
+            const vector& Ci = localCoarseCf[coarseFaceI];
+            vector Ain = Ai/mag(Ai);
+            vector R1i = Ci + (mag(Ai)/wideBy2)*(Ain ^ emptyDir);
+            vector R2i = Ci - (mag(Ai)/wideBy2)*(Ain ^ emptyDir) ;
+
+            const label fromPatchId = compactPatchId[coarseFaceI];
+            patchArea[fromPatchId] += mag(Ai);
+
+            const labelList& visCoarseFaces = visibleFaceFaces[coarseFaceI];
+            forAll(visCoarseFaces, visCoarseFaceI)
+            {
+                F[coarseFaceI].setSize(visCoarseFaces.size());
+                label compactJ = visCoarseFaces[visCoarseFaceI];
+                const vector& Aj = compactCoarseSf[compactJ];
+                const vector& Cj = compactCoarseCf[compactJ];
+
+                const label toPatchId = compactPatchId[compactJ];
+
+                vector Ajn = Aj/mag(Aj);
+                vector R1j = Cj + (mag(Aj)/wideBy2)*(Ajn ^ emptyDir);
+                vector R2j = Cj - (mag(Aj)/wideBy2)*(Ajn ^ emptyDir);
+
+                scalar d1 = mag(R1i - R2j);
+                scalar d2 = mag(R2i - R1j);
+                scalar s1 = mag(R1i - R1j);
+                scalar s2 = mag(R2i - R2j);
+
+                scalar Fij = mag((d1 + d2) - (s1 + s2))/(4.0*mag(Ai)/wideBy2);
+
+                F[coarseFaceI][visCoarseFaceI] = Fij;
+                sumViewFactorPatch[fromPatchId][toPatchId] += Fij*mag(Ai);
+            }
+        }
+    }
+
+    if (Pstream::master())
+    {
+        Info << "Writing view factor matrix..." << endl;
    }

    // Write view factors matrix in listlist form
-    F2LI.write();
+    F.write();

    reduce(sumViewFactorPatch, sumOp<scalarSquareMatrix>());
    reduce(patchArea, sumOp<scalarList>());

+
    if (Pstream::master() && debug)
    {
        forAll(viewFactorsPatches, i)
        {
            label patchI =  viewFactorsPatches[i];
-            for (label j=i; j<viewFactorsPatches.size(); j++)
+            forAll(viewFactorsPatches, i)
            {
-                label patchJ =  viewFactorsPatches[j];
-
+                label patchJ =  viewFactorsPatches[i];
                Info << "F" << patchI << patchJ << ": "
                     << sumViewFactorPatch[patchI][patchJ]/patchArea[patchI]
                     << endl;
@ -1120,13 +758,9 @@ int main(int argc, char *argv[])
        }
    }

+
    if (writeViewFactors)
    {
-        if (Pstream::master())
-        {
-            Info << "Writing view factor matrix..." << endl;
-        }
-
        volScalarField viewFactorField
        (
            IOobject
@ -1157,14 +791,13 @@ int main(int argc, char *argv[])

                forAll(coarseToFine, coarseI)
                {
-                    const scalar FiSum = sum(F2LI[compactI]);
-
+                    const scalar Fij = sum(F[compactI]);
                    const label coarseFaceID = coarsePatchFace[coarseI];
                    const labelList& fineFaces = coarseToFine[coarseFaceID];
                    forAll(fineFaces, fineId)
                    {
                        const label faceID = fineFaces[fineId];
-                        vfbf[patchID][faceID] = FiSum;
+                        vfbf[patchID][faceID] = Fij;
                    }
                    compactI++;
                }
--- a/applications/utilities/surface/surfaceBooleanFeatures/Allwmake
+++ b/applications/utilities/surface/surfaceBooleanFeatures/Allwmake
@ -10,7 +10,7 @@ if have_cgal
 then
    wmake $targetType PolyhedronReader
    export COMP_FLAGS='-IPolyhedronReader'
-    export LINK_FLAGS='-lPolyhedronReader'
+    export LINK_FLAGS='${CGAL_LIBS} -lPolyhedronReader'
 else
    export COMP_FLAGS="-DNO_CGAL"
 fi
--- a/applications/utilities/surface/surfaceBooleanFeatures/Make/options
+++ b/applications/utilities/surface/surfaceBooleanFeatures/Make/options
@ -26,5 +26,4 @@ EXE_LIBS = \
    -lfileFormats \
    -lsurfMesh \
    -lmeshTools \
-    ${CGAL_LIBS} \
    $(LINK_FLAGS)
--- a/applications/utilities/surface/surfaceMeshExtract/surfaceMeshExtract.C
+++ b/applications/utilities/surface/surfaceMeshExtract/surfaceMeshExtract.C
@ -120,6 +120,7 @@ int main(int argc, char *argv[])

    // Less frequently used - reduce some clutter
    argList::setAdvanced("decomposeParDict");
+    argList::setAdvanced("noFunctionObjects");

    argList::addArgument("output", "The output surface file");

--- a/applications/utilities/surface/surfaceToPatch/surfaceToPatch.C
+++ b/applications/utilities/surface/surfaceToPatch/surfaceToPatch.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2020-2022 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -58,7 +58,7 @@ label addPatch(polyMesh& mesh, const word& patchName)
    {
        const polyBoundaryMesh& patches = mesh.boundaryMesh();

-        polyPatchList newPatches(patches.size() + 1);
+        List<polyPatch*> newPatches(patches.size() + 1);

        patchi = 0;

@ -67,25 +67,20 @@ label addPatch(polyMesh& mesh, const word& patchName)
        {
            const polyPatch& pp = patches[i];

-            newPatches.set
-            (
-                patchi,
+            newPatches[patchi] =
                pp.clone
                (
                    patches,
                    patchi,
                    pp.size(),
                    pp.start()
-                )
-            );
+                ).ptr();

            patchi++;
        }

        // Add zero-sized patch
-        newPatches.set
-        (
-            patchi,
+        newPatches[patchi] =
            new polyPatch
            (
                patchName,
@ -94,8 +89,7 @@ label addPatch(polyMesh& mesh, const word& patchName)
                patchi,
                patches,
                polyPatch::typeName
-            )
-        );
+            );

        mesh.removeBoundary();
        mesh.addPatches(newPatches);
--- a/applications/utilities/surface/surfaceTransformPoints/surfaceTransformPoints.C
+++ b/applications/utilities/surface/surfaceTransformPoints/surfaceTransformPoints.C
@ -188,18 +188,15 @@ int main(int argc, char *argv[])
    );
    argList::addBoolOption
    (
-        "auto-centre",
-        "Use bounding box centre as centre for rotations"
+        "auto-origin",
+        "Use bounding box centre as origin for rotations"
    );
    argList::addOption
    (
-        "centre",
+        "origin",
        "point",
-        "Use specified <point> as centre for rotations"
+        "Use specified <point> as origin for rotations"
    );
-    argList::addOptionCompat("auto-centre", {"auto-origin", 2206});
-    argList::addOptionCompat("centre", {"origin", 2206});
-
    argList::addOption
    (
        "rotate",
@ -355,18 +352,18 @@ int main(int argc, char *argv[])
        points += v;
    }

-    vector rotationCentre;
-    bool useRotationCentre = args.readIfPresent("centre", rotationCentre);
-    if (args.found("auto-centre") && !useRotationCentre)
+    vector origin;
+    bool useOrigin = args.readIfPresent("origin", origin);
+    if (args.found("auto-origin") && !useOrigin)
    {
-        useRotationCentre = true;
-        rotationCentre = boundBox(points).centre();
+        useOrigin = true;
+        origin = boundBox(points).centre();
    }

-    if (useRotationCentre)
+    if (useOrigin)
    {
-        Info<< "Set centre of rotation to " << rotationCentre << endl;
-        points -= rotationCentre;
+        Info<< "Set origin for rotations to " << origin << endl;
+        points -= origin;
    }


@ -458,15 +455,15 @@ int main(int argc, char *argv[])
        transform(points, rot, points);
    }

-    if (useRotationCentre)
-    {
-        Info<< "Unset centre of rotation from " << rotationCentre << endl;
-        points += rotationCentre;
-    }
-
    // Output scaling
    applyScaling(points, getScalingOpt("write-scale", args));

+    if (useOrigin)
+    {
+        Info<< "Unset origin for rotations from " << origin << endl;
+        points += origin;
+    }
+
    surf1.movePoints(points);
    surf1.write(exportName, writeFileType);

--- a/bin/foamCleanPath
+++ b/bin/foamCleanPath
@ -7,7 +7,7 @@
 #    \\/     M anipulation  |
 #------------------------------------------------------------------------------
 #     Copyright (C) 2011-2016 OpenFOAM Foundation
-#     Copyright (C) 2017-2022 OpenCFD Ltd.
+#     Copyright (C) 2017-2021 OpenCFD Ltd.
 #------------------------------------------------------------------------------
 # License
 #     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@ -48,29 +48,26 @@
 #     - Similarly for c-shell
 #       eval `foamCleanPath -csh-path dir1:dir2`
 #
-# For library paths, it is suggested to use -sh-lib, -env=-lib etc.
-#
-# On Darwin it uses FOAM_LD_LIBRARY_PATH instead of LD_LIBRARY_PATH.
-# This should actually be DYLD_LIBRARY_PATH on Darwin, but setting that
-# or LD_LIBRARY_PATH via a shell-script is disallowed when SIP is active.
+# For library paths, it is suggested to use -sh-lib, or -csh-lib.
+# The will use DYLD_LIBRARY_PATH instead of LD_LIBRARY_PATH on Darwin.
 #
 #------------------------------------------------------------------------------
 printHelp() {
    cat<<USAGE

-Usage: foamCleanPath [OPTION] ENVNAME [filter] .. [filter]
+Usage: foamCleanPath [OPTION] ENV [filter] .. [filter]
       foamCleanPath [OPTION] -env=name [filter] .. [filter]
 options:
  -env=NAME         Evaluate NAME to obtain initial content,
                    Accepts "-env=-path", "-env=-lib" shortcuts for PATH
-                    and LD_LIBRARY_PATH (FOAM_LD_LIBRARY_PATH on Darwin)
+                    and LD_LIBRARY_PATH (DYLD_LIBRARY_PATH on Darwin)
  -sh=NAME          Produce 'NAME=...' output for sh eval
  -csh=NAME         Produce 'setenv NAME ...' output for csh eval
  -sh-env=NAME      Same as -sh=NAME -env=NAME
  -csh-env=NAME     Same as -csh=NAME -env=NAME
  -sh-path | -csh-path  Same as -[c]sh-env=PATH
  -sh-lib  | -csh-lib   Same as -[c]sh-env=LD_LIBRARY_PATH
-                        (FOAM_LD_LIBRARY_PATH on Darwin)
+                        or DYLD_LIBRARY_PATH on Darwin
  -debug            Print debug information to stderr
  -strip            Remove inaccessible directories
  -verbose          Report some progress (input, output, ...)
@ -84,7 +81,7 @@ Prints its argument (which should be a ':' separated list) cleansed from
 Exit status
    0  on success
    1  general error
-    2  initial value of ENVNAME is empty
+    2  initial value of 'path' is empty

 USAGE
    exit 0  # A clean exit
@ -106,69 +103,63 @@ die()
 #-------------------------------------------------------------------------------

 # Input and outputs
-unset dirList shellFlavour shellOutput
+unset dirList shellOutput shellFlavour
 unset optDebug optEnvName optStrip optVerbose

 # Parse options
 while [ "$#" -gt 0 ]
 do
    case "$1" in
-    (-h | -help*)
+    -h | -help*)
        printHelp
        ;;

-    (-csh-lib | -csh-path | -sh-lib | -sh-path)
+    -env=*)
+        name="${1#*=}"
+        [ -n "$name" ] || die "Option '$1' missing an ENVNAME"
+
+        # Accept aliases
+        case "$name" in
+        (-path) name='PATH';;
+        (-lib)
+            case "$(uname -s 2>/dev/null)" in
+                (Darwin) name='DYLD_LIBRARY_PATH';;
+                (*) name='LD_LIBRARY_PATH';;
+            esac
+            ;;
+         esac
+         optEnvName="$name"     # Use for input evaluation
+         ;;
+
+    -csh-path | -sh-path)
        shellFlavour="$1"
-        case "$1" in
-        (*-lib)
-            name='LD_LIBRARY_PATH'
-            if [ "$(uname -s 2>/dev/null)" = Darwin ]
-            then
-                name='FOAM_LD_LIBRARY_PATH' # Shadow DYLD_LIBRARY_PATH
-            fi
-            ;;
-        (*-path)
-            name='PATH'
-            ;;
+        name='PATH'
+        optEnvName="$name"      # Use for input evaluation
+        shellOutput="$name"     # Use for output
+        ;;
+
+    -csh-lib | -sh-lib)
+        shellFlavour="$1"
+        case "$(uname -s 2>/dev/null)" in
+            (Darwin) name='DYLD_LIBRARY_PATH';;
+            (*) name='LD_LIBRARY_PATH';;
        esac
        optEnvName="$name"      # Use for input evaluation
        shellOutput="$name"     # Use for output
        ;;

-    (-env=*)
-        name="${1#*=}"
-        [ -n "$name" ] || die "Option '$1' missing an ENVNAME"
-        # Handle (-lib | -path) aliases
-        case "$1" in
-        (*=-lib)
-            name='LD_LIBRARY_PATH'
-            if [ "$(uname -s 2>/dev/null)" = Darwin ]
-            then
-                name='FOAM_LD_LIBRARY_PATH' # Shadow DYLD_LIBRARY_PATH
-            fi
-            ;;
-        (*=-path)
-            name='PATH'
-            ;;
-        esac
-        optEnvName="$name"      # Use for input evaluation
-        ;;
-
-    (-csh=* | -csh-env=* | -sh=* | -sh-env=*)
+    -csh=* | -sh=* | -csh-env=* | -sh-env=*)
        shellFlavour="$1"
        name="${1#*=}"
        [ -n "$name" ] || die "Option '$1' missing an ENVNAME"
-        # Handle (-lib | -path) aliases
-        case "$1" in
-        (*=-lib)
-            name='LD_LIBRARY_PATH'
-            if [ "$(uname -s 2>/dev/null)" = Darwin ]
-            then
-                name='FOAM_LD_LIBRARY_PATH' # Shadow DYLD_LIBRARY_PATH
-            fi
-            ;;
-        (*=-path)
-            name='PATH'
+        # Accept aliases
+        case "$name" in
+        (-path) name='PATH';;
+        (-lib)
+            case "$(uname -s 2>/dev/null)" in
+                (Darwin) name='DYLD_LIBRARY_PATH';;
+                (*) name='LD_LIBRARY_PATH';;
+            esac
            ;;
        esac
        shellOutput="$name"     # Use for output
@ -176,16 +167,16 @@ do
        case "$1" in (*-env=*) optEnvName="$name";; esac
        ;;

-    (-debug)
+    -debug)
        optDebug=true
        ;;
-    (-strip)
+    -strip)
        optStrip=true
        ;;
-    (-verbose)
+    -verbose)
        optVerbose=true
        ;;
-    (*)
+    *)
        break
        ;;
    esac
--- a/Show More
+++ b/Show More