BUG: integratedNonUniformTable: correct offsets (fixes #2614)

- with geometryOrder=1, calculate the edge normals from the adjacent
  faces (area-weighted, inverse distance squared) and also
  use that for the Le() calculation.

  Includes the contributions from processor edge neighbours, so it
  should be consistent on both sides.

  This new method (consider as 'beta') contrasts with the current
  standard method that first calculates area-weighted point normals
  and uses the average of them for the edge normal.

  Enable for testing either with a controlDict OptimisationSwitch entry
  "fa:geometryOrder", or on the command-line:

      solverName -opt-switch=fa:geometryOrder=1

META-INFO/api-info

@@ -1,2 +1,2 @@
 api=2206
-patch=0
+patch=220907

applications/solvers/acoustic/acousticFoam/Make/options

@@ -11,4 +11,5 @@ EXE_LIBS = \
     -lfvOptions \
     -lmeshTools \
     -lsampling \
-    -lregionFaModels
+    -lregionFaModels \
+    -lfiniteArea

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 \

applications/solvers/compressible/rhoPimpleFoam/Make/options

@@ -26,4 +26,5 @@ EXE_LIBS = \
     -ltopoChangerFvMesh \
     -lsampling \
     -latmosphericModels \
-    -lregionFaModels
+    -lregionFaModels \
+    -lfiniteArea

applications/solvers/compressible/rhoSimpleFoam/Make/options

@@ -21,4 +21,5 @@ EXE_LIBS = \
     -lcompressibleTurbulenceModels \
     -lthermoTools \
     -latmosphericModels \
-    -lregionFaModels
+    -lregionFaModels \
+    -lfiniteArea

applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options

@@ -28,4 +28,5 @@ EXE_LIBS = \
     -ldynamicMesh \
     -ldynamicFvMesh \
     -ltopoChangerFvMesh \
-    -lregionFaModels
+    -lregionFaModels \
+    -lfiniteArea

applications/solvers/heatTransfer/buoyantSimpleFoam/Make/options

@@ -22,4 +22,5 @@ EXE_LIBS = \
     -lcompressibleTurbulenceModels \
     -lthermoTools \
     -latmosphericModels \
-    -lregionFaModels
+    -lregionFaModels \
+    -lfiniteArea

applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options

@@ -40,4 +40,5 @@ EXE_LIBS = \
     -lthermoTools \
     -lradiationModels \
     -lregionModels \
-    -lregionFaModels
+    -lregionFaModels \
+    -lfiniteArea

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/Make/options

@@ -36,5 +36,5 @@ EXE_LIBS = \
     -lfaOptions \
     -lregionModels \
     -lsampling \
-    -lregionFaModels
-
+    -lregionFaModels \
+    -lfiniteArea

applications/solvers/incompressible/pimpleFoam/Make/options

@@ -22,4 +22,5 @@ EXE_LIBS = \
     -ldynamicFvMesh \
     -ltopoChangerFvMesh \
     -latmosphericModels \
-    -lregionFaModels
+    -lregionFaModels \
+    -lfiniteArea

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 \

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 \

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 \

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 \

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 \

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 \

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 \

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 \

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 \

applications/solvers/multiphase/MPPICInterFoam/Make/options

@@ -42,4 +42,5 @@ EXE_LIBS = \
     -lregionModels \
     -lregionFaModels \
     -lsurfaceFilmModels \
+    -lfiniteArea \
     -lcompressibleTwoPhaseMixtureTurbulenceModels

applications/test/DynamicList2/Test-DynamicList2.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2021-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -104,7 +104,7 @@ void readList
 {
     OTstream os;
     os << input;
-    ITstream is("input", os.tokens());
+    ITstream is(os.tokens());
 
     is >> output;
 }
@@ -118,7 +118,7 @@ void readList
 {
     OTstream os;
     os << input;
-    ITstream is("input", os.tokens());
+    ITstream is(os.tokens());
 
     is >> output;
 }

applications/test/FieldFields1/Make/files

@@ -0,0 +1,3 @@
+Test-FieldFields1.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-FieldFields1

applications/test/FieldFields1/Test-FieldFields1.C

@@ -0,0 +1,117 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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;
+}
+
+
+// ************************************************************************* //

applications/test/FieldFields2/Make/files

@@ -0,0 +1,3 @@
+Test-FieldFields2.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-FieldFields2

applications/test/FieldFields2/Make/options

@@ -0,0 +1,2 @@
+/* EXE_INC = */
+/* EXE_LIBS = */

applications/test/FieldFields2/Test-FieldFields2.C

@@ -24,7 +24,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    Test-tensorFieldFields1
+    Test-FieldFields2
 
 \*---------------------------------------------------------------------------*/
 
@@ -73,21 +73,29 @@ void allocComponents
 }
 
 
-vectorField randomVectorField(label size)
+template<class Type>
+tmp<Field<Type>> randomField(Random& rnd, label dim)
 {
-    Random rnd;
+    auto tfld = tmp<Field<Type>>::New(dim);
+    auto& fld = tfld.ref();
 
-    vectorField vf(size);
-
-    forAll(vf, i)
+    for (Type& val : fld)
     {
-        for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt)
+        for (direction cmpt=0; cmpt < pTraits<Type>::nComponents; ++cmpt)
         {
-            vf[i][cmpt] = rnd.position<label>(0, 100);
+            setComponent(val, cmpt) = rnd.position<label>(0, 100);
         }
     }
 
-    return vf;
+    return tfld;
+}
+
+
+template<class Type>
+tmp<Field<Type>> randomField(label dim)
+{
+    Random rnd;
+    return randomField<Type>(rnd, dim);
 }
 
 
@@ -191,7 +199,7 @@ int main(int argc, char *argv[])
         printFieldField(sf1);
 
         Info<< nl;
-        for (direction cmpt = 0; cmpt < vector::nComponents; ++cmpt)
+        for (direction cmpt = 0; cmpt < pTraits<vector>::nComponents; ++cmpt)
         {
             unzipRow(sf1, cmpt, slice[0]);
 
@@ -253,8 +261,9 @@ int main(int argc, char *argv[])
     {
         Info<< nl << "vectorField" << nl;
 
+        Random rnd;
         FieldField<Field, vector> vf1(1);
-        vf1.set(0, new vectorField(randomVectorField(4)));
+        vf1.set(0, randomField<vector>(rnd, 4));
 
         FixedList<FieldField<Field, scalar>, 3> cmpts;
         allocComponents(cmpts, 4);

applications/test/HashPtrTable/Test-HashPtrTable.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,7 +31,10 @@ Description
 #include <memory>
 #include <iostream>
 #include "autoPtr.H"
-#include "HashPtrTable.H"
+#include "refPtr.H"
+#include "tmp.H"
+#include "PtrMap.H"
+#include "primitiveFields.H"
 
 using namespace Foam;
 
@@ -250,6 +253,42 @@ 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;
 }
 

applications/test/Hashing1/Make/options

@@ -1,3 +1,3 @@
-EXE_INC = ${c++LESSWARN}
+EXE_INC = $(c++LESSWARN)
 
 /* EXE_LIBS = */

applications/test/Hashing2/Make/options

@@ -1,3 +1,3 @@
-EXE_INC = ${c++LESSWARN}
+EXE_INC = $(c++LESSWARN)
 
 /* EXE_LIBS = */

applications/test/HashingSpeed/Make/options

@@ -1,3 +1,3 @@
-EXE_INC = ${c++LESSWARN}
+EXE_INC = $(c++LESSWARN)
 
 /* EXE_LIBS = */

applications/test/Polynomial/Test-Polynomial.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2022 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("input", std::move(os.tokens()));
+        ITstream is(std::move(os.tokens()));
         is >> polyfunc;
     }
 

applications/test/TestTools/TestTools.H

@@ -28,15 +28,14 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-using namespace Foam;
-
-#include "floatScalar.H"
-#include "doubleScalar.H"
+#include "scalar.H"
 #include "complex.H"
 #include "Matrix.H"
 #include "Random.H"
 #include <chrono>
 
+using namespace Foam;
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 // Total number of unit tests

applications/test/autoPtr/Test-autoPtr.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,7 +27,6 @@ License
 
 // #define Foam_autoPtr_deprecate_setMethod
 
-#include <memory>
 #include "autoPtr.H"
 #include "labelList.H"
 #include "ListOps.H"
@@ -77,6 +76,14 @@ autoPtr<T> testNullReturn2()
 }
 
 
+template<class T>
+struct DerivedList : public List<T>
+{
+    // Inherit constructors
+    using List<T>::List;
+};
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 //  Main program:
 
@@ -254,6 +261,16 @@ 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;
 }

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;
 
-    defineTemplate2TypeNameAndDebug(dimUint8Field, 0);
+    defineTemplateTypeNameAndDebug(dimUint8Field, 0);
 
 } // End namespace Foam
 #endif

applications/test/dynamicLibrary/Test-dynamicLibrary.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020-2021 OpenCFD Ltd.
+    Copyright (C) 2020-2022 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");
 

applications/test/faceHashing/Make/options

@@ -1,3 +1,3 @@
-EXE_INC = ${c++LESSWARN}
+EXE_INC = $(c++LESSWARN)
 
 /* EXE_LIBS = */

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 \

applications/test/fvc2D/Test-fvc2D.C

@@ -39,7 +39,7 @@ namespace Foam
 {
     typedef GeometricField<vector2D, fvPatchField, volMesh> volVector2DField;
 
-    defineTemplate2TypeNameAndDebug(volVector2DField::Internal, 0);
+    defineTemplateTypeNameAndDebug(volVector2DField::Internal, 0);
     defineTemplateTypeNameAndDebug(volVector2DField, 0);
 
     typedef fvPatchField<vector2D> fvPatchVector2DField;

applications/test/matrices/DiagonalMatrix/Test-DiagonalMatrix.C

@@ -38,8 +38,7 @@ Description
 
 #include "DiagonalMatrix.H"
 #include "RectangularMatrix.H"
-#include "floatScalar.H"
-#include "doubleScalar.H"
+#include "scalar.H"
 #include "complex.H"
 #include "TestTools.H"
 

applications/test/matrices/RectangularMatrix/Test-RectangularMatrix.C

@@ -44,8 +44,7 @@ Note
 
 #include "RectangularMatrix.H"
 #include "SquareMatrix.H"
-#include "floatScalar.H"
-#include "doubleScalar.H"
+#include "scalar.H"
 #include "complex.H"
 #include "IOmanip.H"
 #include "TestTools.H"

applications/test/matrices/SquareMatrix/Test-SquareMatrix.C

@@ -45,8 +45,7 @@ Note
 #include "scalarMatrices.H"
 #include "RectangularMatrix.H"
 #include "SquareMatrix.H"
-#include "floatScalar.H"
-#include "doubleScalar.H"
+#include "scalar.H"
 #include "complex.H"
 #include "IOmanip.H"
 #include "TestTools.H"

applications/test/matrices/SymmetricSquareMatrix/Test-SymmetricSquareMatrix.C

@@ -46,8 +46,7 @@ Note
 #include "RectangularMatrix.H"
 #include "SquareMatrix.H"
 #include "SymmetricSquareMatrix.H"
-#include "floatScalar.H"
-#include "doubleScalar.H"
+#include "scalar.H"
 #include "complex.H"
 #include "IOmanip.H"
 #include "Random.H"

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) */

applications/test/parallel-chunks/Test-parallel-chunks.C

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

applications/test/refPtr/Test-refPtr.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020-2021 OpenCFD Ltd.
+    Copyright (C) 2020-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -19,6 +19,9 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "primitiveFields.H"
+#include "autoPtr.H"
+#include "refPtr.H"
+#include "tmp.H"
 #include "Switch.H"
 
 using namespace Foam;
@@ -29,6 +32,36 @@ 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)
 {
@@ -37,15 +70,24 @@ void printInfo(const refPtr<T>& item, const bool verbose = false)
         << " addr: " << Foam::name(item.get())
         << " movable:" << Switch(item.movable());
 
-    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;
+    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>>();
 
     if (verbose && item)
     {
@@ -101,6 +143,62 @@ 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;

applications/test/scalarOps/Test-scalarOps.C

@@ -126,6 +126,9 @@ int main(int argc, char *argv[])
     testDivide<vector>(vectors);
 
 
+    //(void) compareOp<vector>()(vector::zero, vector::one);
+
+
     Info<< "\nEnd\n" << endl;
 
     return 0;

applications/test/tensorFieldFields1/Make/files

@@ -1,3 +0,0 @@
-Test-tensorFieldFields1.C
-
-EXE = $(FOAM_USER_APPBIN)/Test-tensorFieldFields1

applications/test/tensorFields1/Test-tensorFields1.C

@@ -34,21 +34,29 @@ Application
 
 using namespace Foam;
 
-vectorField randomVectorField(label size)
+template<class Type>
+tmp<Field<Type>> randomField(Random& rnd, label dim)
 {
-    Random rnd;
+    auto tfld = tmp<Field<Type>>::New(dim);
+    auto& fld = tfld.ref();
 
-    vectorField vf(size);
-
-    forAll(vf, i)
+    for (Type& val : fld)
     {
-        for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt)
+        for (direction cmpt=0; cmpt < pTraits<Type>::nComponents; ++cmpt)
         {
-            vf[i][cmpt] = rnd.position<label>(0, 100);
+            setComponent(val, cmpt) = rnd.position<label>(0, 100);
         }
     }
 
-    return vf;
+    return tfld;
+}
+
+
+template<class Type>
+tmp<Field<Type>> randomField(label dim)
+{
+    Random rnd;
+    return randomField<Type>(rnd, dim);
 }
 
 
@@ -251,7 +259,7 @@ int main(int argc, char *argv[])
     {
         Info<< nl << "vectorField" << nl;
 
-        vectorField vf1(randomVectorField(4));
+        vectorField vf1(randomField<vector>(4));
         FixedList<scalarField, 3> cmpts(scalarField(vf1.size()));
 
         Info<< nl

applications/test/timeSelector/Test-timeSelector.C

@@ -63,7 +63,6 @@ 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");

applications/test/tmp/Test-tmp.C

@@ -19,6 +19,9 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "primitiveFields.H"
+#include "autoPtr.H"
+#include "refPtr.H"
+#include "tmp.H"
 #include "Switch.H"
 
 using namespace Foam;
@@ -126,6 +129,28 @@ 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;

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;
 
-    defineTemplate2TypeNameAndDebug(volUint8Field, 0);
+    defineTemplateTypeNameAndDebug(volUint8Field, 0);
 
 } // End namespace Foam
 #endif

applications/utilities/finiteArea/checkFaMesh/printMeshSummary.H

@@ -116,12 +116,16 @@ Description
     {
         const faPatch& p = patches[patchi];
 
-         Info<< "  " << "patch " << p.index()
-             << " (size: " << returnReduce(p.size(), sumOp<label>())
-             << ") name: " << p.name()
-             << nl;
+        // Report physical size (nEdges) not virtual size
+        Info<< "  " << "patch " << p.index()
+            << " (size: " << returnReduce(p.nEdges(), sumOp<label>())
+            << ") name: " << p.name()
+            << nl;
     }
 
+    Info<< "----------------" << nl
+        << "Used polyPatches: " << flatOutput(aMesh.whichPolyPatches()) << nl;
+
 
     // Geometry information
     Info<< nl;

applications/utilities/finiteArea/makeFaMesh/printMeshSummary.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2021-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -36,12 +36,17 @@ Description
     {
         const faPatch& p = patches[patchi];
 
-         Info<< "  " << "patch " << p.index()
-             << " (size: " << returnReduce(p.size(), sumOp<label>())
-             << ") name: " << p.name()
-             << nl;
+        // Report physical size (nEdges) not virtual size
+        Info<< "  " << "patch " << p.index()
+            << " (size: " << returnReduce(p.nEdges(), sumOp<label>())
+            << ") name: " << p.name()
+            << nl;
     }
 
+    Info<< "----------------" << nl
+        << "Used polyPatches: " << flatOutput(aMesh.whichPolyPatches()) << nl;
+
+
     // Geometry information
     Info<< nl;
     {

applications/utilities/mesh/advanced/refinementLevel/refinementLevel.C

@@ -226,14 +226,18 @@ int main(int argc, char *argv[])
     // Add the boundary patches
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
-    List<polyPatch*> p(patches.size());
+    polyPatchList newPatches(patches.size());
 
-    forAll(p, patchi)
+    forAll(newPatches, patchi)
     {
-        p[patchi] = patches[patchi].clone(fMesh.boundaryMesh()).ptr();
+        newPatches.set
+        (
+            patchi,
+            patches[patchi].clone(fMesh.boundaryMesh())
+        );
     }
 
-    fMesh.addFvPatches(p);
+    fMesh.addFvPatches(newPatches);
 
 
     // Refinement level

applications/utilities/mesh/advanced/snappyRefineMesh/snappyRefineMesh.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2019 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -205,34 +205,43 @@ label addPatch(polyMesh& mesh, const word& patchName)
     {
         const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
-        List<polyPatch*> newPatches(patches.size() + 1);
+        polyPatchList newPatches(patches.size() + 1);
+
+        label nPatches = 0;
 
         // Add empty patch as 0th entry (Note: only since subsetMesh wants this)
         patchi = 0;
-
-        newPatches[patchi] =
-            new emptyPolyPatch
-            (
-                Foam::word(patchName),
-                0,
-                mesh.nInternalFaces(),
-                patchi,
-                patches,
-                emptyPolyPatch::typeName
-            );
-
-        forAll(patches, i)
         {
-            const polyPatch& pp = patches[i];
+            newPatches.set
+            (
+                nPatches,
+                new emptyPolyPatch
+                (
+                    patchName,
+                    0,
+                    mesh.nInternalFaces(),
+                    nPatches,
+                    patches,
+                    emptyPolyPatch::typeName
+                )
+            );
+            ++nPatches;
+        }
 
-            newPatches[i+1] =
+        for (const polyPatch& pp : patches)
+        {
+            newPatches.set
+            (
+                nPatches,
                 pp.clone
                 (
                     patches,
-                    i+1,
+                    nPatches,
                     pp.size(),
                     pp.start()
-                ).ptr();
+                )
+            );
+            ++nPatches;
         }
 
         mesh.removeBoundary();

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)
     {
-        List<polyPatch*> newPatchPtrList((mesh.boundaryMesh().size() - 1));
+        polyPatchList newPatches((mesh.boundaryMesh().size() - 1));
         label newPatchi = 0;
         forAll(mesh.boundaryMesh(), patchi)
         {
@@ -1686,18 +1686,21 @@ int main(int argc, char *argv[])
             {
                 const polyPatch& patch = mesh.boundaryMesh()[patchi];
 
-                newPatchPtrList[newPatchi] = patch.clone
+                newPatches.set
                 (
-                    mesh.boundaryMesh(),
                     newPatchi,
-                    patch.size(),
-                    patch.start()
-                ).ptr();
-
-                newPatchi++;
+                    patch.clone
+                    (
+                        mesh.boundaryMesh(),
+                        newPatchi,
+                        patch.size(),
+                        patch.start()
+                    )
+                );
+                ++newPatchi;
             }
         }
-        repatcher.changePatches(newPatchPtrList);
+        repatcher.changePatches(newPatches);
     }
 
     // Set the precision of the points data to 10

applications/utilities/mesh/generation/blockMesh/mergePatchPairs.H

@@ -1,30 +1,56 @@
-// Handle merging of patch pairs
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+\*---------------------------------------------------------------------------*/
+
 {
     wordPairList mergePatchPairs;
 
-    // Read in a list of dictionaries for the merge patch pairs
+    // Read in a list of merge patch pairs
     if
     (
         meshDict.readIfPresent("mergePatchPairs", mergePatchPairs)
      && mergePatchPairs.size()
     )
     {
-        Info<< "Creating merge patch pairs" << nl << endl;
+        Info<< "Merging " << mergePatchPairs.size() << " patch pairs" << nl;
 
-        Info<< "Adding point and face zones" << endl;
+        // Cleanup
+        wordHashSet cleanupPatches(4*mergePatchPairs.size());
+        wordHashSet cleanupPointZones(2*mergePatchPairs.size());
+        wordHashSet cleanupFaceZones(2*mergePatchPairs.size());
+
+        Info<< "    Adding point and face zones" << endl;
         {
-            auto& pzs = mesh.pointZones();
-            pzs.clearAddressing();
-            auto& fzs = mesh.faceZones();
-            fzs.clearAddressing();
+            const auto& pbm = mesh.boundaryMesh();
+
+            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()
-                  + name(pairi)
+                  + Foam::name(pairi)
                 );
 
                 // An empty zone for cut points
@@ -37,40 +63,35 @@
                         pzs
                     )
                 );
+                cleanupPointZones.insert(pzs.last().name());
 
-                // Master patch
-                const word masterPatchName(mergePatchPairs[pairi].first());
-                const polyPatch& masterPatch =
-                    mesh.boundaryMesh()[masterPatchName];
-
+                // Coupling side 0 (master)
                 fzs.append
                 (
                     new faceZone
                     (
-                        mergeName + "MasterZone",
-                        identity(masterPatch.range()),
+                        mergeName + "Side0Zone",
+                        identity(patch0.range()),
                         false, // none are flipped
                         fzs.size(),
                         fzs
                     )
                 );
+                cleanupFaceZones.insert(fzs.last().name());
 
-                // Slave patch
-                const word slavePatchName(mergePatchPairs[pairi].second());
-                const polyPatch& slavePatch =
-                    mesh.boundaryMesh()[slavePatchName];
-
+                // Coupling side 1 (slave)
                 fzs.append
                 (
                     new faceZone
                     (
-                        mergeName + "SlaveZone",
-                        identity(slavePatch.range()),
+                        mergeName + "Side1Zone",
+                        identity(patch1.range()),
                         false, // none are flipped
                         fzs.size(),
                         fzs
                     )
                 );
+                cleanupFaceZones.insert(fzs.last().name());
 
                 // An empty zone for cut faces
                 fzs.append
@@ -82,35 +103,38 @@
                         fzs
                     )
                 );
-            }  // end of all merge pairs
+                cleanupFaceZones.insert(fzs.last().name());
+            }
         }
 
 
-
-        Info<< "Creating attachPolyTopoChanger" << endl;
+        Info<< "    Merging with attachPolyTopoChanger" << endl;
         attachPolyTopoChanger polyMeshAttacher(mesh);
-        polyMeshAttacher.setSize(mergePatchPairs.size());
+        polyMeshAttacher.resize(1);
 
         forAll(mergePatchPairs, pairi)
         {
+            cleanupPatches.insert(mergePatchPairs[pairi].first());
+            cleanupPatches.insert(mergePatchPairs[pairi].second());
+
             const word mergeName
             (
                 mergePatchPairs[pairi].first()
               + mergePatchPairs[pairi].second()
-              + name(pairi)
+              + Foam::name(pairi)
             );
 
             // Add the sliding interface mesh modifier
             polyMeshAttacher.set
             (
-                pairi,
+                0,
                 new slidingInterface
                 (
-                    "couple" + name(pairi),
+                    "couple" + Foam::name(pairi),
                     pairi,
                     polyMeshAttacher,
-                    mergeName + "MasterZone",
-                    mergeName + "SlaveZone",
+                    mergeName + "Side0Zone",
+                    mergeName + "Side1Zone",
                     mergeName + "CutPointZone",
                     mergeName + "CutFaceZone",
                     mergePatchPairs[pairi].first(),
@@ -120,12 +144,135 @@
                     intersection::VISIBLE
                 )
             );
+
+            polyMeshAttacher.attach(false);  // Do not yet remove empty patches
         }
 
-        polyMeshAttacher.attach(true);
+        // 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;
+            }
+        }
     }
     else
     {
-        Info<< nl << "There are no merge patch pairs" << endl;
+        Info<< "No patch pairs to merge" << endl;
     }
 }
+
+
+// ************************************************************************* //

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)
         {

applications/utilities/mesh/generation/foamyMesh/foamyQuadMesh/foamyQuadMesh.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2021-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -167,28 +167,32 @@ int main(int argc, char *argv[])
     );
 
     Info<< "Constructing patches." << endl;
-    List<polyPatch*> patches(poly2DMesh.patchNames().size());
-    label countPatches = 0;
+    polyPatchList newPatches(poly2DMesh.patchNames().size());
+    label nPatches = 0;
 
-    forAll(patches, patchi)
+    forAll(newPatches, patchi)
     {
         if (poly2DMesh.patchSizes()[patchi] != 0)
         {
-            patches[countPatches] = new polyPatch
+            newPatches.set
             (
-                poly2DMesh.patchNames()[patchi],
-                poly2DMesh.patchSizes()[patchi],
-                poly2DMesh.patchStarts()[patchi],
-                countPatches,
-                pMesh.boundaryMesh(),
-                word::null
+                nPatches,
+                new polyPatch
+                (
+                    poly2DMesh.patchNames()[patchi],
+                    poly2DMesh.patchSizes()[patchi],
+                    poly2DMesh.patchStarts()[patchi],
+                    nPatches,
+                    pMesh.boundaryMesh(),
+                    word::null
+                )
             );
 
-            countPatches++;
+            ++nPatches;
         }
     }
-    patches.setSize(countPatches);
-    pMesh.addPatches(patches);
+    newPatches.resize(nPatches);
+    pMesh.addPatches(newPatches);
 
     if (extrude)
     {

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C

@@ -1342,6 +1342,17 @@ 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
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

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
-    List<polyPatch*> newPatchPtrList(patches.size());
+    polyPatchList newPatches(patches.size());
 
     newPatchi = 0;
 
@@ -200,34 +200,41 @@ int main(int argc, char *argv[])
     {
         const polyPatch& patch = mesh.boundaryMesh()[patchi];
 
-        newPatchPtrList[newPatchi] =
+        newPatches.set
+        (
+            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 boundaryPatch& bp = patches[patchi];
+        const word& patchName = patches[patchi].name();
 
-        newPatchPtrList[newPatchi] = polyPatch::New
+        newPatches.set
         (
-            polyPatch::typeName,
-            bp.name(),
-            0,
-            mesh.nFaces(),
             newPatchi,
-            mesh.boundaryMesh()
-        ).ptr();
+            polyPatch::New
+            (
+                polyPatch::typeName,
+                patchName,
+                0,
+                mesh.nFaces(),
+                newPatchi,
+                mesh.boundaryMesh()
+            )
+        );
 
-        newPatchi++;
+        ++newPatchi;
     }
 
     if (!overwrite)
@@ -238,7 +245,7 @@ int main(int argc, char *argv[])
 
     // Change patches
     repatchPolyTopoChanger polyMeshRepatcher(mesh);
-    polyMeshRepatcher.changePatches(newPatchPtrList);
+    polyMeshRepatcher.changePatches(newPatches);
 
 
     // Change face ordering

applications/utilities/mesh/manipulation/checkMesh/checkTopology.C

@@ -40,6 +40,7 @@ License
 #include "vtkSurfaceWriter.H"
 #include "checkTools.H"
 #include "treeBoundBox.H"
+#include "syncTools.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -48,25 +49,113 @@ 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) << returnReduce(pp.size(), sumOp<label>())
+        << setw(9) << globalSize
         << setw(9) << returnReduce(pp.nPoints(), sumOp<label>());
 
-    if (!Pstream::parRun())
+    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
     {
-        typedef typename PatchType::surfaceTopo TopoType;
         TopoType pTyp = pp.surfaceType();
 
-        if (pp.empty())
-        {
-            Info<< setw(34) << "ok (empty)";
-        }
-        else if (pTyp == TopoType::MANIFOLD)
+        if (pTyp == TopoType::MANIFOLD)
         {
             if (pp.checkPointManifold(true, &points))
             {
@@ -549,15 +638,8 @@ Foam::label Foam::checkTopology
     );
 
     {
-        if (!Pstream::parRun())
-        {
-            Info<< "\nChecking patch topology for multiply connected"
+        Info<< "\nChecking patch topology for multiply connected"
                 << " surfaces..." << endl;
-        }
-        else
-        {
-            Info<< "\nChecking basic patch addressing..." << endl;
-        }
 
         const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
@@ -566,11 +648,8 @@ Foam::label Foam::checkTopology
         Info<< "    "
             << setw(20) << "Patch"
             << setw(9) << "Faces"
-            << setw(9) << "Points";
-        if (!Pstream::parRun())
-        {
-            Info<< setw(34) << "Surface topology";
-        }
+            << setw(9) << "Points"
+            << "Surface topology";
         if (allGeometry)
         {
             Info<< " Bounding box";
@@ -583,7 +662,16 @@ Foam::label Foam::checkTopology
 
             if (!isA<processorPolyPatch>(pp))
             {
-                checkPatch(allGeometry, pp.name(), pp, points);
+                checkPatch
+                (
+                    allGeometry,
+                    pp.name(),
+                    mesh,
+                    pp,
+                    identity(pp.size(), pp.start()),
+                    pp.meshEdges(),
+                    points
+                );
                 Info<< endl;
             }
         }
@@ -592,15 +680,8 @@ Foam::label Foam::checkTopology
     }
 
     {
-        if (!Pstream::parRun())
-        {
-            Info<< "\nChecking faceZone topology for multiply connected"
-                << " surfaces..." << endl;
-        }
-        else
-        {
-            Info<< "\nChecking basic faceZone addressing..." << endl;
-        }
+        Info<< "\nChecking faceZone topology for multiply connected"
+            << " surfaces..." << endl;
 
         Pout.setf(ios_base::left);
 
@@ -611,12 +692,8 @@ Foam::label Foam::checkTopology
             Info<< "    "
                 << setw(20) << "FaceZone"
                 << setw(9) << "Faces"
-                << setw(9) << "Points";
-
-            if (!Pstream::parRun())
-            {
-                Info<< setw(34) << "Surface topology";
-            }
+                << setw(9) << "Points"
+                << setw(34) << "Surface topology";
             if (allGeometry)
             {
                 Info<< " Bounding box";
@@ -625,7 +702,16 @@ Foam::label Foam::checkTopology
 
             for (const faceZone& fz : faceZones)
             {
-                checkPatch(allGeometry, fz.name(), fz(), points);
+                checkPatch
+                (
+                    allGeometry,
+                    fz.name(),
+                    mesh,
+                    fz(),           // patch
+                    fz,             // mesh face labels
+                    fz.meshEdges(), // mesh edge labels
+                    points
+                );
                 Info<< endl;
             }
         }

applications/utilities/mesh/manipulation/checkMesh/checkTopology.H

@@ -1,4 +1,4 @@
-#include "label.H"
+#include "labelList.H"
 #include "autoPtr.H"
 
 namespace Foam
@@ -14,7 +14,10 @@ namespace Foam
     (
         const bool allGeometry,
         const word& name,
+        const polyMesh& mesh,
         const PatchType& pp,
+        const labelList& meshFaces,
+        const labelList& meshEdges,
         pointSet& points
     );
 

applications/utilities/mesh/manipulation/mergeMeshes/mergePolyMesh.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2016 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -136,46 +136,55 @@ void Foam::mergePolyMesh::sortProcessorPatches()
 
     const polyBoundaryMesh& oldPatches = boundaryMesh();
 
-    DynamicList<polyPatch*> newPatches(oldPatches.size());
+    polyPatchList newPatches(oldPatches.size());
 
     labelList oldToSorted(oldPatches.size());
 
+    label nPatches = 0;
+
     forAll(oldPatches, patchi)
     {
         const polyPatch& pp = oldPatches[patchi];
 
         if (!isA<processorPolyPatch>(pp))
         {
-            oldToSorted[patchi] = newPatches.size();
-            newPatches.append
+            newPatches.set
             (
+                nPatches,
                 pp.clone
                 (
                     oldPatches,
-                    oldToSorted[patchi],
+                    nPatches,
                     0,
                     nInternalFaces()
-                ).ptr()
+                )
             );
+
+            oldToSorted[patchi] = nPatches;
+            ++nPatches;
         }
     }
+
     forAll(oldPatches, patchi)
     {
         const polyPatch& pp = oldPatches[patchi];
 
         if (isA<processorPolyPatch>(pp))
         {
-            oldToSorted[patchi] = newPatches.size();
-            newPatches.append
+            newPatches.set
             (
+                nPatches,
                 pp.clone
                 (
                     oldPatches,
                     oldToSorted[patchi],
                     0,
                     nInternalFaces()
-                ).ptr()
+                )
             );
+
+            oldToSorted[patchi] = nPatches;
+            ++nPatches;
         }
     }
 
@@ -461,7 +470,7 @@ void Foam::mergePolyMesh::merge()
     {
         Info<< "Copying old patches" << endl;
 
-        List<polyPatch*> newPatches(patchNames_.size());
+        polyPatchList newPatches(patchNames_.size());
 
         const polyBoundaryMesh& oldPatches = boundaryMesh();
 
@@ -470,7 +479,11 @@ void Foam::mergePolyMesh::merge()
 
         for (patchi = 0; patchi < oldPatches.size(); patchi++)
         {
-            newPatches[patchi] = oldPatches[patchi].clone(oldPatches).ptr();
+            newPatches.set
+            (
+                patchi,
+                oldPatches[patchi].clone(oldPatches)
+            );
         }
 
         Info<< "Adding new patches. " << endl;
@@ -487,15 +500,18 @@ void Foam::mergePolyMesh::merge()
             dict.set("nFaces", 0);
             dict.set("startFace", endOfLastPatch);
 
-            newPatches[patchi] =
+            newPatches.set
             (
-                polyPatch::New
+                patchi,
                 (
-                    patchNames_[patchi],
-                    dict,
-                    patchi,
-                    oldPatches
-                ).ptr()
+                    polyPatch::New
+                    (
+                        patchNames_[patchi],
+                        dict,
+                        patchi,
+                        oldPatches
+                    )
+                )
             );
         }
 

applications/utilities/mesh/manipulation/mirrorMesh/mirrorFvMesh.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2017-2018 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -404,20 +404,24 @@ Foam::mirrorFvMesh::mirrorFvMesh
     fvMesh& pMesh = mirrorMeshPtr_();
 
     // Add the boundary patches
-    List<polyPatch*> p(newPatchSizes.size());
+    polyPatchList newPatches(newPatchSizes.size());
 
-    forAll(p, patchi)
+    forAll(newPatches, patchi)
     {
-        p[patchi] = boundaryMesh()[patchi].clone
+        newPatches.set
         (
-            pMesh.boundaryMesh(),
             patchi,
-            newPatchSizes[patchi],
-            newPatchStarts[patchi]
-        ).ptr();
+            boundaryMesh()[patchi].clone
+            (
+                pMesh.boundaryMesh(),
+                patchi,
+                newPatchSizes[patchi],
+                newPatchStarts[patchi]
+            )
+        );
     }
 
-    pMesh.addPatches(p);
+    pMesh.addPatches(newPatches);
 }
 
 

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 \

applications/utilities/mesh/manipulation/stitchMesh/stitchMesh.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2017-2020 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -564,7 +564,7 @@ int main(int argc, char *argv[])
 
             mesh.faceZones()
             (
-                mergePatchName + "MasterZone",
+                mergePatchName + "Side0Zone",
                 true // verbose
             ).resetAddressing(std::move(faceIds), false);
 
@@ -574,7 +574,7 @@ int main(int argc, char *argv[])
 
             mesh.faceZones()
             (
-                mergePatchName + "SlaveZone",
+                mergePatchName + "Side1Zone",
                 true // verbose
             ).resetAddressing(std::move(faceIds), false);
 
@@ -595,8 +595,8 @@ int main(int argc, char *argv[])
                     "couple" + Foam::name(actioni),
                     0,
                     stitcher,
-                    mergePatchName + "MasterZone",
-                    mergePatchName + "SlaveZone",
+                    mergePatchName + "Side0Zone",
+                    mergePatchName + "Side1Zone",
                     mergePatchName + "CutPointZone",
                     cutZoneName,
                     masterPatchName,

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

@@ -62,8 +62,8 @@ Usage
     -rotate-z angle
         Rotate (degrees) about z-axis.
 
-     or -yawPitchRoll (yawdegrees pitchdegrees rolldegrees)
-     or -rollPitchYaw (rolldegrees pitchdegrees yawdegrees)
+     or -yawPitchRoll : (yaw pitch roll) degrees
+     or -rollPitchYaw : (roll pitch yaw) degrees
 
     -scale scalar|vector
         Scale the points by the given scalar or vector on output.
@@ -268,15 +268,18 @@ int main(int argc, char *argv[])
     );
     argList::addBoolOption
     (
-        "auto-origin",
-        "Use bounding box centre as origin for rotations"
+        "auto-centre",
+        "Use bounding box centre as centre for rotations"
     );
     argList::addOption
     (
-        "origin",
+        "centre",
         "point",
-        "Use specified <point> as origin for rotations"
+        "Use specified <point> as centre for rotations"
     );
+    argList::addOptionCompat("auto-centre", {"auto-origin", 2206});
+    argList::addOptionCompat("centre", {"origin", 2206});
+
     argList::addOption
     (
         "rotate",
@@ -437,18 +440,18 @@ int main(int argc, char *argv[])
             points += v;
         }
 
-        vector origin;
-        bool useOrigin = args.readIfPresent("origin", origin);
-        if (args.found("auto-origin") && !useOrigin)
+        vector rotationCentre;
+        bool useRotationCentre = args.readIfPresent("centre", rotationCentre);
+        if (args.found("auto-centre") && !useRotationCentre)
         {
-            useOrigin = true;
-            origin = boundBox(points).centre();
+            useRotationCentre = true;
+            rotationCentre = boundBox(points).centre();
         }
 
-        if (useOrigin)
+        if (useRotationCentre)
         {
-            Info<< "Set origin for rotations to " << origin << endl;
-            points -= origin;
+            Info<< "Set centre of rotation to " << rotationCentre << endl;
+            points -= rotationCentre;
         }
 
 
@@ -545,15 +548,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()));

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 pOpen(const std::string&, int line=0);
+static std::string pipeOpen(const std::string& cmd, const int lineNum = 0);
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -132,51 +132,50 @@ void version()
 }
 
 
-std::string pOpen(const std::string& cmd, int line)
+// 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 res;
+    std::string str;
 
-    FILE* cmdPipe = popen(cmd.c_str(), "r");
-    if (!cmdPipe) return res;
+    FILE* handle = popen(cmd.c_str(), "r");
+    if (!handle) return str;
 
     char* buf = nullptr;
+    size_t len = 0;
+    ssize_t nread;
 
-    // Read line number of lines
-    for (int cnt = 0; cnt <= line; ++cnt)
+    // Read lineNum number of lines
+    for
+    (
+        int cnt = 0;
+        cnt <= lineNum && (nread = ::getline(&buf, &len, handle)) >= 0;
+        ++cnt
+    )
     {
-        size_t linecap = 0;
-        ssize_t linelen = ::getline(&buf, &linecap, cmdPipe);
-
-        if (linelen < 0)
+        if (cnt == lineNum)
         {
-            break;
-        }
+            // Retain the last line, trimming trailing newline
+            str.assign(buf);
 
-        if (cnt == line)
-        {
-            res = std::string(buf);
-
-            // Trim trailing newline
-            if (res.size())
+            if (str.size())
             {
-                res.resize(res.size()-1);
+                str.resize(str.size()-1);
             }
-            break;
         }
     }
 
-    if (buf) free(buf);
+    free(buf);
+    pclose(handle);
 
-    pclose(cmdPipe);
-
-    return res;
+    return str;
 }
 
 
 std::string getLine(const std::string& filename, const std::string& addr)
 {
     std::string line =
-        pOpen
+        pipeOpen
         (
             "echo 'image lookup -va " + addr
           + "'"

applications/utilities/miscellaneous/foamHasLibrary/foamHasLibrary.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020-2021 OpenCFD Ltd.
+    Copyright (C) 2020-2022 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",

applications/utilities/parallelProcessing/decomposePar/decomposePar.C

@@ -1014,7 +1014,9 @@ int main(int argc, char *argv[])
                                 Time::controlDictName,
                                 args.rootPath(),
                                 args.caseName()
-                              / ("processor" + Foam::name(proci))
+                              / ("processor" + Foam::name(proci)),
+                                args.allowFunctionObjects(),
+                                args.allowLibs()
                             )
                         );
                     }

applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C

@@ -244,7 +244,9 @@ int main(int argc, char *argv[])
             (
                 Time::controlDictName,
                 args.rootPath(),
-                args.caseName()/("processor" + Foam::name(proci))
+                args.caseName()/("processor" + Foam::name(proci)),
+                args.allowFunctionObjects(),
+                args.allowLibs()
             )
         );
     }

applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C

@@ -731,7 +731,9 @@ int main(int argc, char *argv[])
             (
                 Time::controlDictName,
                 args.rootPath(),
-                args.caseName()/("processor" + Foam::name(proci))
+                args.caseName()/("processor" + Foam::name(proci)),
+                args.allowFunctionObjects(),
+                args.allowLibs()
             )
         );
     }
@@ -1184,7 +1186,7 @@ int main(int argc, char *argv[])
                         false
                     );
 
-                    masterMeshPtr = fvMeshes[0];
+                    masterMeshPtr.cref(fvMeshes[0]);
                 }
 
 

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

@@ -168,7 +168,6 @@ int main(int argc, char *argv[])
 
     // Less frequently used - reduce some clutter
     argList::setAdvanced("decomposeParDict");
-    argList::setAdvanced("noFunctionObjects");
 
     argList::addVerboseOption("Additional verbosity");
 

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

@@ -65,7 +65,7 @@ label writeDimFields
 {
     typedef VolumeInternalField<Type> FieldType;
 
-    const fvMesh& mesh = dynamicCast<const fvMesh>(ensMesh.mesh());
+    const auto& mesh = refCast<const fvMesh>(ensMesh.mesh());
 
     label count = 0;
 

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

@@ -88,7 +88,7 @@ label writeVolFields
 {
     typedef VolumeField<Type> FieldType;
 
-    const fvMesh& mesh = dynamicCast<const fvMesh>(ensMesh.mesh());
+    const auto& mesh = refCast<const fvMesh>(ensMesh.mesh());
 
     label count = 0;
 

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

@@ -232,13 +232,13 @@ vtk::outputOptions getOutputOptions(const argList& args)
 
         if (!args.found("ascii"))
         {
-            if (sizeof(floatScalar) != 4 || sizeof(label) != 4)
+            if (sizeof(float) != 4 || sizeof(label) != 4)
             {
                 opts.ascii(true);
 
                 WarningInFunction
                     << "Using ASCII rather than legacy binary VTK format since "
-                    << "floatScalar and/or label are not 4 bytes in size."
+                    << "float and/or label are not 4 bytes in size."
                     << nl << endl;
             }
             else
@@ -268,7 +268,6 @@ int main(int argc, char *argv[])
 
     // Less frequently used - reduce some clutter
     argList::setAdvanced("decomposeParDict");
-    argList::setAdvanced("noFunctionObjects");
 
     argList::addVerboseOption("Additional verbosity");
 

applications/utilities/postProcessing/lagrangian/particleTracks/particleTracks.C

@@ -87,7 +87,6 @@ int main(int argc, char *argv[])
 
     // Less frequently used - reduce some clutter
     argList::setAdvanced("decomposeParDict");
-    argList::setAdvanced("noFunctionObjects");
 
     argList::addOption
     (

applications/utilities/postProcessing/miscellaneous/profilingSummary/profilingSummary.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2020 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -100,7 +100,9 @@ int main(int argc, char *argv[])
             (
                 Time::controlDictName,
                 args.rootPath(),
-                args.caseName()/("processor" + Foam::name(proci))
+                args.caseName()/("processor" + Foam::name(proci)),
+                args.allowFunctionObjects(),
+                args.allowLibs()
             )
         );
     }

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 = IOdictionary::typeName;
-                const_cast<word&>(IOdictionary::typeName) = word::null;
+                const word oldTypeName = localIOdictionary::typeName;
+                const_cast<word&>(localIOdictionary::typeName) = word::null;
 
                 IOobject fieldHeader
                 (
@@ -675,11 +675,15 @@ int main(int argc, char *argv[])
                     false
                 );
 
-                if (fieldHeader.typeHeaderOk<IOdictionary>(false))
+                if (fieldHeader.typeHeaderOk<localIOdictionary>(false))
                 {
-                    IOdictionary fieldDict(fieldHeader);
+                    //IOdictionary fieldDict(fieldHeader);
+                    //- dictionaries to-be-changed are either boundary
+                    //  or field dictionary. Both are processor-local.
+                    localIOdictionary fieldDict(fieldHeader);
 
-                    const_cast<word&>(IOdictionary::typeName) = oldTypeName;
+                    const_cast<word&>(localIOdictionary::typeName) =
+                        oldTypeName;
 
                     // Fake type back to what was in field
                     const_cast<word&>(fieldDict.type()) =

applications/utilities/preProcessing/faceAgglomerate/faceAgglomerate.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2020,2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -96,7 +96,15 @@ int main(int argc, char *argv[])
 
     label nCoarseFaces = 0;
 
-    for (const entry& dEntry : agglomDict)
+
+    const auto& patchesDict =
+        agglomDict.optionalSubDict
+        (
+            "patchAgglomeration",
+            keyType::LITERAL
+        );
+
+    for (const entry& dEntry : patchesDict)
     {
         labelList patchids = boundary.indices(dEntry.keyword());
 
@@ -112,7 +120,7 @@ int main(int argc, char *argv[])
                 (
                     pp.localFaces(),
                     pp.localPoints(),
-                    agglomDict.subDict(pp.name())
+                    dEntry.dict()
                 );
 
                 agglomObject.agglomerate();

applications/utilities/preProcessing/mapFieldsPar/MapVolFields.H

@@ -43,7 +43,8 @@ void evaluateConstraintTypes(GeometricField<Type, fvPatchField, volMesh>& fld)
 {
     auto& fldBf = fld.boundaryFieldRef();
 
-    const UPstream::commsTypes commsType(UPstream::defaultCommsType);
+    const UPstream::commsTypes commsType = UPstream::defaultCommsType;
+    const label startOfRequests = UPstream::nRequests();
 
     if
     (
@@ -51,8 +52,6 @@ void evaluateConstraintTypes(GeometricField<Type, fvPatchField, volMesh>& fld)
      || commsType == UPstream::commsTypes::nonBlocking
     )
     {
-        const label startOfRequests = UPstream::nRequests();
-
         forAll(fldBf, patchi)
         {
             fvPatchField<Type>& tgtField = fldBf[patchi];
@@ -68,11 +67,7 @@ void evaluateConstraintTypes(GeometricField<Type, fvPatchField, volMesh>& fld)
         }
 
         // Wait for outstanding requests
-        if
-        (
-            UPstream::parRun()
-         && commsType == UPstream::commsTypes::nonBlocking
-        )
+        if (commsType == UPstream::commsTypes::nonBlocking)
         {
             UPstream::waitRequests(startOfRequests);
         }

applications/utilities/preProcessing/setExprBoundaryFields/setExprBoundaryFields.C

@@ -56,6 +56,7 @@ 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/

applications/utilities/preProcessing/setExprFields/setExprFields.C

@@ -494,6 +494,7 @@ 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/
@@ -741,7 +742,7 @@ int main(int argc, char *argv[])
             );
         }
 
-        if (args.found("withFunctionObjects"))
+        if (args.allowFunctionObjects())
         {
             runTime.functionObjects().start();
         }

applications/utilities/preProcessing/setFields/Make/options

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

applications/utilities/preProcessing/setFields/setFields.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,214 +38,311 @@ 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 word& fieldTypeDesc,
+    const fieldDescription& fieldDesc,
     const fvMesh& mesh,
     const labelList& selectedCells,
-    Istream& fieldValueStream
+    Istream& is
 )
 {
     typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
 
-    if (fieldTypeDesc != fieldType::typeName + "Value")
+    if (fieldDesc.type() != fieldType::typeName)
     {
         return false;
     }
 
-    word fieldName(fieldValueStream);
+    // Get value from stream
+    const Type fieldValue = pTraits<Type>(is);
+
 
     // Check the current time directory
     IOobject fieldHeader
     (
-        fieldName,
-        mesh.time().timeName(),
-        mesh,
+        fieldDesc.name(),
+        mesh.thisDb().time().timeName(),
+        mesh.thisDb(),
         IOobject::MUST_READ
     );
 
-    // Check the "constant" directory
-    if (!fieldHeader.typeHeaderOk<fieldType>(true))
+    bool found = fieldHeader.typeHeaderOk<fieldType>(true);
+
+    if (!found)
     {
+        // Fallback to "constant" directory
         fieldHeader = IOobject
         (
-            fieldName,
-            mesh.time().constant(),
-            mesh,
+            fieldDesc.name(),
+            mesh.thisDb().time().constant(),
+            mesh.thisDb(),
             IOobject::MUST_READ
         );
+        found = fieldHeader.typeHeaderOk<fieldType>(true);
     }
 
-    // Check field exists
-    if (fieldHeader.typeHeaderOk<fieldType>(true))
+    // Field exists
+    if (found)
     {
-        Info<< "    Setting internal values of "
+        Info<< "    - set internal values of "
             << fieldHeader.headerClassName()
-            << " " << fieldName << endl;
+            << ": " << fieldDesc.name()
+            << " = " << fieldValue << endl;
 
         fieldType field(fieldHeader, mesh, false);
 
-        const Type value = pTraits<Type>(fieldValueStream);
-
-        if (selectedCells.size() == field.size())
+        if (isNull(selectedCells) || selectedCells.size() == field.size())
         {
-            field.primitiveFieldRef() = value;
+            field.primitiveFieldRef() = fieldValue;
         }
         else
         {
-            forAll(selectedCells, celli)
+            for (const label celli : selectedCells)
             {
-                field[selectedCells[celli]] = value;
+                field[celli] = fieldValue;
             }
         }
 
-        typename GeometricField<Type, fvPatchField, volMesh>::
-            Boundary& fieldBf = field.boundaryFieldRef();
-
-        forAll(field.boundaryField(), patchi)
+        // Make boundary fields consistent - treat like zeroGradient
+        for (auto& pfld : field.boundaryFieldRef())
         {
-            fieldBf[patchi] = fieldBf[patchi].patchInternalField();
+            pfld = pfld.patchInternalField();
         }
 
         if (!field.write())
         {
             FatalErrorInFunction
-              << "Failed writing field " << fieldName << endl;
+                << "Failed writing field " << field.name() << endl;
         }
     }
     else
     {
-        WarningInFunction
-          << "Field " << fieldName << " not found" << endl;
-
-        // Consume value
-        (void)pTraits<Type>(fieldValueStream);
+        Warning
+            << "Field " << fieldDesc.name() << " not found" << endl;
     }
 
     return true;
 }
 
 
-class setCellField
-{
-
-public:
-
-    setCellField()
-    {}
-
-    autoPtr<setCellField> clone() const
-    {
-        return autoPtr<setCellField>::New();
-    }
-
-    class iNew
-    {
-        const fvMesh& mesh_;
-        const labelList& selectedCells_;
-
-    public:
-
-        iNew(const fvMesh& mesh, const labelList& selectedCells)
-        :
-            mesh_(mesh),
-            selectedCells_(selectedCells)
-        {}
-
-        iNew(const fvMesh& mesh, labelList&& selectedCells)
-        :
-            mesh_(mesh),
-            selectedCells_(std::move(selectedCells))
-        {}
-
-        autoPtr<setCellField> operator()(Istream& fieldValues) const
-        {
-            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)
-                )
-            )
-            {
-                WarningInFunction
-                    << "field type " << fieldType << " not currently supported"
-                    << endl;
-            }
-
-            return autoPtr<setCellField>::New();
-        }
-    };
-};
-
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// Setting finite-area face fields
 template<class Type>
-bool setFaceFieldType
+bool setAreaFieldType
 (
-    const word& fieldTypeDesc,
-    const fvMesh& mesh,
+    const fieldDescription& fieldDesc,
+    const faMesh& mesh,
     const labelList& selectedFaces,
-    Istream& fieldValueStream
+    Istream& is
 )
 {
-    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+    typedef GeometricField<Type, faPatchField, areaMesh> fieldType;
 
-    if (fieldTypeDesc != fieldType::typeName + "Value")
+    if (fieldDesc.type() != fieldType::typeName)
     {
         return false;
     }
 
-    word fieldName(fieldValueStream);
+    // Get value from stream
+    const Type fieldValue = pTraits<Type>(is);
+
 
     // Check the current time directory
     IOobject fieldHeader
     (
-        fieldName,
-        mesh.time().timeName(),
-        mesh,
+        fieldDesc.name(),
+        mesh.thisDb().time().timeName(),
+        mesh.thisDb(),
         IOobject::MUST_READ
     );
 
-    // Check the "constant" directory
-    if (!fieldHeader.typeHeaderOk<fieldType>(true))
+    bool found = fieldHeader.typeHeaderOk<fieldType>(true);
+
+    if (!found)
     {
+        // Fallback to "constant" directory
         fieldHeader = IOobject
         (
-            fieldName,
-            mesh.time().constant(),
-            mesh,
+            fieldDesc.name(),
+            mesh.thisDb().time().constant(),
+            mesh.thisDb(),
             IOobject::MUST_READ
         );
+        found = fieldHeader.typeHeaderOk<fieldType>(true);
     }
 
-    // Check field exists
-    if (fieldHeader.typeHeaderOk<fieldType>(true))
+    // Field exists
+    if (found)
     {
-        Info<< "    Setting patchField values of "
+        Info<< "    - set internal values of "
             << fieldHeader.headerClassName()
-            << " " << fieldName << endl;
+            << ": " << fieldDesc.name()
+            << " = " << fieldValue << endl;
 
         fieldType field(fieldHeader, mesh);
 
-        const Type value = pTraits<Type>(fieldValueStream);
+        if (isNull(selectedFaces) || selectedFaces.size() == field.size())
+        {
+            field.primitiveFieldRef() = fieldValue;
+        }
+        else
+        {
+            for (const label facei : selectedFaces)
+            {
+                field[facei] = fieldValue;
+            }
+        }
+
+        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 fvMesh& mesh,
+    const labelList& selectedFaces,
+    Istream& is
+)
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+
+    if (fieldDesc.type() != fieldType::typeName)
+    {
+        return false;
+    }
+
+    // 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)
+    {
+        // Fallback to "constant" directory
+        fieldHeader = IOobject
+        (
+            fieldDesc.name(),
+            mesh.thisDb().time().constant(),
+            mesh.thisDb(),
+            IOobject::MUST_READ
+        );
+        found = fieldHeader.typeHeaderOk<fieldType>(true);
+    }
+
+    // Field exists
+    if (found)
+    {
+        Info<< "    - set boundary values of "
+            << fieldHeader.headerClassName()
+            << ": " << fieldDesc.name()
+            << " = " << fieldValue << endl;
+
+        fieldType field(fieldHeader, mesh);
 
         // Create flat list of selected faces and their value.
         Field<Type> allBoundaryValues(mesh.nBoundaryFaces());
@@ -260,34 +358,47 @@ bool setFaceFieldType
         }
 
         // Override
-        bool hasWarned = false;
+        unsigned hasWarned = 0;
         labelList nChanged
         (
             returnReduce(field.boundaryField().size(), maxOp<label>()),
-            0
+            Zero
         );
-        forAll(selectedFaces, i)
+
+        for (const label facei : selectedFaces)
         {
-            label facei = selectedFaces[i];
-            if (mesh.isInternalFace(facei))
+            const label bFacei = facei-mesh.nInternalFaces();
+
+            if (bFacei < 0)
             {
-                if (!hasWarned)
+                if (!(hasWarned & 1))
                 {
-                    hasWarned = true;
+                    hasWarned |= 1;
                     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 bFacei = facei-mesh.nInternalFaces();
-                allBoundaryValues[bFacei] = value;
-                nChanged[mesh.boundaryMesh().patchID()[bFacei]]++;
+                label patchi = mesh.boundaryMesh().patchID()[bFacei];
+
+                allBoundaryValues[bFacei] = fieldValue;
+                ++nChanged[patchi];
             }
         }
 
-        Pstream::listCombineAllGather(nChanged, plusEqOp<label>());
+        Pstream::listCombineReduce(nChanged, plusEqOp<label>());
 
         auto& fieldBf = field.boundaryFieldRef();
 
@@ -299,6 +410,7 @@ bool setFaceFieldType
                 Info<< "    On patch "
                     << field.boundaryField()[patchi].patch().name()
                     << " set " << nChanged[patchi] << " values" << endl;
+
                 fieldBf[patchi] == SubField<Type>
                 (
                     allBoundaryValues,
@@ -312,85 +424,222 @@ bool setFaceFieldType
         if (!field.write())
         {
             FatalErrorInFunction
-                << "Failed writing field " << field.name() << exit(FatalError);
+                << "Failed writing field " << field.name() << endl;
         }
     }
     else
     {
-        WarningInFunction
-          << "Field " << fieldName << " not found" << endl;
-
-        // Consume value
-        (void)pTraits<Type>(fieldValueStream);
+        Warning
+            << "Field " << fieldDesc.name() << " not found" << endl;
     }
 
     return true;
 }
 
 
-class setFaceField
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Dispatcher for setting volume fields
+struct setCellField
 {
+    autoPtr<setCellField> clone() const { return nullptr; }  // placeholder
 
-public:
-
-    setFaceField()
-    {}
-
-    autoPtr<setFaceField> clone() const
+    static bool apply
+    (
+        const fieldDescription& fieldDesc,
+        const fvMesh& m,
+        const labelList& selectedCells,
+        Istream& is
+    )
     {
-        return autoPtr<setFaceField>::New();
+        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)
+        );
     }
 
     class iNew
     {
         const fvMesh& mesh_;
-        const labelList& selectedFaces_;
+        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_;
 
     public:
 
         iNew(const fvMesh& mesh, const labelList& selectedFaces)
         :
             mesh_(mesh),
-            selectedFaces_(selectedFaces)
+            selected_(selectedFaces)
         {}
 
-        iNew(const fvMesh& mesh, labelList&& selectedFaces)
-        :
-            mesh_(mesh),
-            selectedFaces_(std::move(selectedFaces))
-        {}
-
-        autoPtr<setFaceField> operator()(Istream& fieldValues) const
+        autoPtr<setFaceField> operator()(Istream& is) const
         {
-            word fieldType(fieldValues);
+            const fieldDescription fieldDesc(is);
 
-            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)
-                )
-            )
+            bool ok = setFaceField::apply(fieldDesc, mesh_, selected_, is);
+
+            if (!ok)
             {
-                WarningInFunction
-                    << "field type " << fieldType << " not currently supported"
-                    << endl;
+                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 autoPtr<setFaceField>::New();
+            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)
+        :
+            mesh_(mesh),
+            selected_(selectedFaces)
+        {}
+
+        autoPtr<setAreaField> operator()(Istream& is) const
+        {
+            const fieldDescription fieldDesc(is);
+
+            bool ok = setAreaField::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
+        }
+    };
+};
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -403,11 +652,28 @@ 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"
 
@@ -415,70 +681,144 @@ int main(int argc, char *argv[])
 
     IOdictionary setFieldsDict(dictIO);
 
-    if (setFieldsDict.found("defaultFieldValues"))
+
+    // 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
     {
-        Info<< "Setting field default values" << endl;
-        PtrList<setCellField> defaultFieldValues
-        (
-            setFieldsDict.lookup("defaultFieldValues"),
-            setCellField::iNew(mesh, labelList(mesh.nCells()))
-        );
-        Info<< endl;
+        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 region values" << endl;
+    Info<< "Setting field region values" << nl << endl;
 
     PtrList<entry> regions(setFieldsDict.lookup("regions"));
 
-    forAll(regions, regionI)
+    for (const entry& region : regions)
     {
-        const entry& region = regions[regionI];
-
         autoPtr<topoSetSource> source =
             topoSetSource::New(region.keyword(), mesh, region.dict());
 
         if (source().setType() == topoSetSource::CELLSET_SOURCE)
         {
-            cellSet selectedCellSet
+            cellSet subset
             (
                 mesh,
                 "cellSet",
                 mesh.nCells()/10+1  // Reasonable size estimate.
             );
 
-            source->applyToSet
-            (
-                topoSetSource::NEW,
-                selectedCellSet
-            );
+            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");
 
             PtrList<setCellField> fieldValues
             (
-                region.dict().lookup("fieldValues"),
-                setCellField::iNew(mesh, selectedCellSet.sortedToc())
+                is,
+                setCellField::iNew(mesh, selectedCells)
             );
         }
         else if (source().setType() == topoSetSource::FACESET_SOURCE)
         {
-            faceSet selectedFaceSet
+            faceSet subset
             (
                 mesh,
                 "faceSet",
                 mesh.nBoundaryFaces()/10+1
             );
 
-            source->applyToSet
-            (
-                topoSetSource::NEW,
-                selectedFaceSet
-            );
+            source->applyToSet(topoSetSource::NEW, subset);
+
+            labelList selectedFaces(subset.sortedToc());
+
+            Info<< "    Selected " << selectedFaces.size()
+                << " faces" << nl;
+
+            ITstream& is = region.dict().lookup("fieldValues");
 
             PtrList<setFaceField> fieldValues
             (
-                region.dict().lookup("fieldValues"),
-                setFaceField::iNew(mesh, selectedFaceSet.sortedToc())
+                is,
+                setFaceField::iNew(mesh, selectedFaces)
             );
+
+            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)
+                );
+            }
         }
     }
 

applications/utilities/preProcessing/setTurbulenceFields/setTurbulenceFields.C

@@ -109,6 +109,7 @@ Note
 #include "singlePhaseTransportModel.H"
 #include "turbulentTransportModel.H"
 #include "turbulentFluidThermoModel.H"
+#include "processorFvPatchField.H"
 #include "wallFvPatch.H"
 #include "fixedValueFvPatchFields.H"
 
@@ -122,6 +123,40 @@ 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,
@@ -447,22 +482,26 @@ 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())
@@ -475,6 +514,7 @@ int main(int argc, char *argv[])
         {
             R[celli] = Rdiag[celli];
         }
+        correctProcessorPatches<symmTensor>(R);
     }
 
 

applications/utilities/preProcessing/viewFactorsGen/Allwmake

@@ -0,0 +1,19 @@
+#!/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
+
+#------------------------------------------------------------------------------

applications/utilities/preProcessing/viewFactorsGen/Make/options

@@ -1,10 +1,18 @@
+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 \

applications/utilities/preProcessing/viewFactorsGen/searchingEngine_CGAL.H

@@ -0,0 +1,84 @@
+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());

applications/utilities/preProcessing/viewFactorsGen/shootRays.H

@@ -1,6 +1,15 @@
 // 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
 (

applications/utilities/preProcessing/viewFactorsGen/shootRays_CGAL.H

@@ -0,0 +1,131 @@
+// 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();
+}

applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C

@@ -31,14 +31,34 @@ Group
     grpPreProcessingUtilities
 
 Description
-    View factors are calculated based on a face agglomeration array
-    (finalAgglom generated by faceAgglomerate utility).
+    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.
 
-    Each view factor between the agglomerated faces i and j (Fij) is calculated
-    using a double integral of the sub-areas composing the agglomeration.
+    The input from viewFactorsDict are:
 
-    The patches involved in the view factor calculation are taken from the
-    boundary file and should be part on the group viewFactorWall. ie.:
+        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.
 
     floor
     {
@@ -48,6 +68,9 @@ 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"
@@ -56,21 +79,52 @@ 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"
 
-using namespace Foam;
+#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
 (
@@ -127,20 +181,56 @@ triSurface triangulate
         newPatchI++;
     }
 
-    //striSurfaceToAgglom.resize(localTriFaceI-1);
+    //triSurfaceToAgglom.resize(localTriFaceI-1);
 
     triangles.shrink();
+    triSurface surface(triangles, mesh.points());
+    surface.compactPoints();
 
-    // Create globally numbered tri surface
-    triSurface rawSurface(triangles, mesh.points());
 
-    // Create locally numbered tri surface
-    triSurface surface
+#ifndef NO_CGAL
+
+    // CGAL : every processor has whole surface
+
+    globalIndex globalFaceIdx(surface.size(), globalIndex::gatherOnly());
+    globalIndex globalPointIdx
     (
-        rawSurface.localFaces(),
-        rawSurface.localPoints()
+        surface.points().size(), globalIndex::gatherOnly()
     );
 
+    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);
 
@@ -193,7 +283,7 @@ void writeRays
 }
 
 
-scalar calculateViewFactorFij
+scalar calculateViewFactorFij2AI
 (
     const vector& i,
     const vector& j,
@@ -250,6 +340,74 @@ 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[])
@@ -273,7 +431,7 @@ int main(int argc, char *argv[])
             "viewFactorsDict",
             runTime.constant(),
             mesh,
-            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::MUST_READ,
             IOobject::NO_WRITE
        )
     );
@@ -288,6 +446,23 @@ 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
     (
@@ -296,12 +471,22 @@ int main(int argc, char *argv[])
             "finalAgglom",
             mesh.facesInstance(),
             mesh,
-            IOobject::MUST_READ,
+            IOobject::READ_IF_PRESENT,
             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
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -336,10 +521,8 @@ 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();
@@ -465,7 +648,13 @@ 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
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -473,11 +662,15 @@ 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);
@@ -518,9 +711,9 @@ int main(int argc, char *argv[])
         visibleFaceFaces[faceI][nVisibleFaceFaces[faceI]++] = compactI;
     }
 
-
     // Construct data in compact addressing
-    // I need coarse Sf (Ai), fine Sf (dAi) and fine Cf(r) to calculate Fij
+    // (2AA) need coarse (Ai), fine Sf (dAi) and fine Cf(r) to calculate Fij
+    // (2LI) need edges (li)
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     pointField compactCoarseCf(map.constructSize(), Zero);
@@ -528,12 +721,16 @@ 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)
@@ -548,11 +745,21 @@ 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++];
+                List<point>& fineSf = compactFineSf[compactI];
+
+                label startFace = pp.start();
+
+                const vectorField locPoints
+                (
+                    mesh.points(),
+                    faces[coarseI + startFace]
+                );
 
                 const label coarseFaceI = coarsePatchFace[coarseI];
                 const labelList& fineFaces = coarseToFine[coarseFaceI];
@@ -560,6 +767,8 @@ int main(int argc, char *argv[])
                 fineCf.setSize(fineFaces.size());
                 fineSf.setSize(fineFaces.size());
 
+                compactPoints.append(locPoints);
+
                 fineCf = UIndirectList<point>
                 (
                     mesh.Cf().boundaryField()[patchID],
@@ -570,19 +779,25 @@ 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)
     {
@@ -598,8 +813,7 @@ int main(int argc, char *argv[])
 
     // Fill local view factor matrix
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    scalarListIOList F
+    scalarListIOList F2LI
     (
         IOobject
         (
@@ -630,6 +844,61 @@ 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)
@@ -638,119 +907,212 @@ 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)
             {
-                F[coarseFaceI].setSize(visCoarseFaces.size());
+                //F2AI[coarseFaceI].setSize(visCoarseFaces.size());
+                F2LI[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];
 
-                scalar Fij = 0;
+                bool far(false);
+                // Relative distance
                 forAll(localFineSf, i)
                 {
-                    const vector& dAi = localFineSf[i];
+                    const scalar dAi =
+                        Foam::sqrt
+                        (
+                            mag(localFineSf[i])/constant::mathematical::pi
+                        );
                     const vector& dCi = localFineCf[i];
 
                     forAll(remoteFineSj, j)
                     {
-                        const vector& dAj = remoteFineSj[j];
+                        const scalar dAj =
+                            Foam::sqrt
+                            (
+                                mag(remoteFineSj[j])/constant::mathematical::pi
+                            );
                         const vector& dCj = remoteFineCj[j];
 
-                        scalar dIntFij = calculateViewFactorFij
-                        (
-                            dCi,
-                            dCj,
-                            dAi,
-                            dAj
-                        );
+                        const scalar dist = mag(dCi - dCj)/((dAi + dAj)/2);
 
-                        Fij += dIntFij;
+                        if (dist > distTol)
+                        {
+                            far = true;
+                        }
                     }
                 }
-                F[coarseFaceI][visCoarseFaceI] = Fij/mag(Ai);
-                sumViewFactorPatch[fromPatchId][toPatchId] += Fij;
+
+                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
+                        (
+                            mag(F2LIij) > ROOTVSMALL
+                          ? (F2LIij-oldEToeInt)/F2LIij
+                          : 0
+                        );
+
+                        if
+                        (
+                            (mag(err) < GaussQuadTol && gi > 0)
+                         || gi == maxQuadOrder-1
+                        )
+                        {
+                            F2LI[coarseFaceI][visCoarseFaceI] =
+                                F2LIij/mag(Ai)/4/constant::mathematical::pi;
+                            break;
+                        }
+                        else
+                        {
+                            oldEToeInt = F2LIij;
+                        }
+                    }
+                }
+
+                sumViewFactorPatch[fromPatchId][toPatchId] +=
+                    F2LI[coarseFaceI][visCoarseFaceI]*mag(Ai);
             }
         }
     }
-    else if (mesh.nSolutionD() == 2)
+    else
     {
-        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;
+         FatalErrorInFunction
+            << " View factors are not available in 2D "
+            << exit(FatalError);
     }
 
     // Write view factors matrix in listlist form
-    F.write();
+    F2LI.write();
 
     reduce(sumViewFactorPatch, sumOp<scalarSquareMatrix>());
     reduce(patchArea, sumOp<scalarList>());
 
-
     if (Pstream::master() && debug)
     {
         forAll(viewFactorsPatches, i)
         {
             label patchI =  viewFactorsPatches[i];
-            forAll(viewFactorsPatches, i)
+            for (label j=i; j<viewFactorsPatches.size(); j++)
             {
-                label patchJ =  viewFactorsPatches[i];
+                label patchJ =  viewFactorsPatches[j];
+
                 Info << "F" << patchI << patchJ << ": "
                      << sumViewFactorPatch[patchI][patchJ]/patchArea[patchI]
                      << endl;
@@ -758,9 +1120,13 @@ int main(int argc, char *argv[])
         }
     }
 
-
     if (writeViewFactors)
     {
+        if (Pstream::master())
+        {
+            Info << "Writing view factor matrix..." << endl;
+        }
+
         volScalarField viewFactorField
         (
             IOobject
@@ -791,13 +1157,14 @@ int main(int argc, char *argv[])
 
                 forAll(coarseToFine, coarseI)
                 {
-                    const scalar Fij = sum(F[compactI]);
+                    const scalar FiSum = sum(F2LI[compactI]);
+
                     const label coarseFaceID = coarsePatchFace[coarseI];
                     const labelList& fineFaces = coarseToFine[coarseFaceID];
                     forAll(fineFaces, fineId)
                     {
                         const label faceID = fineFaces[fineId];
-                        vfbf[patchID][faceID] = Fij;
+                        vfbf[patchID][faceID] = FiSum;
                     }
                     compactI++;
                 }

applications/utilities/surface/surfaceBooleanFeatures/Allwmake

@@ -10,7 +10,7 @@ if have_cgal
 then
     wmake $targetType PolyhedronReader
     export COMP_FLAGS='-IPolyhedronReader'
-    export LINK_FLAGS='${CGAL_LIBS} -lPolyhedronReader'
+    export LINK_FLAGS='-lPolyhedronReader'
 else
     export COMP_FLAGS="-DNO_CGAL"
 fi

applications/utilities/surface/surfaceBooleanFeatures/Make/options

@@ -26,4 +26,5 @@ EXE_LIBS = \
     -lfileFormats \
     -lsurfMesh \
     -lmeshTools \
+    ${CGAL_LIBS} \
     $(LINK_FLAGS)

applications/utilities/surface/surfaceMeshExtract/surfaceMeshExtract.C

@@ -120,7 +120,6 @@ 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");
 

applications/utilities/surface/surfaceToPatch/surfaceToPatch.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2020-2021 OpenCFD Ltd.
+    Copyright (C) 2020-2022 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();
 
-        List<polyPatch*> newPatches(patches.size() + 1);
+        polyPatchList newPatches(patches.size() + 1);
 
         patchi = 0;
 
@@ -67,20 +67,25 @@ label addPatch(polyMesh& mesh, const word& patchName)
         {
             const polyPatch& pp = patches[i];
 
-            newPatches[patchi] =
+            newPatches.set
+            (
+                patchi,
                 pp.clone
                 (
                     patches,
                     patchi,
                     pp.size(),
                     pp.start()
-                ).ptr();
+                )
+            );
 
             patchi++;
         }
 
         // Add zero-sized patch
-        newPatches[patchi] =
+        newPatches.set
+        (
+            patchi,
             new polyPatch
             (
                 patchName,
@@ -89,7 +94,8 @@ label addPatch(polyMesh& mesh, const word& patchName)
                 patchi,
                 patches,
                 polyPatch::typeName
-            );
+            )
+        );
 
         mesh.removeBoundary();
         mesh.addPatches(newPatches);

applications/utilities/surface/surfaceTransformPoints/surfaceTransformPoints.C

@@ -188,15 +188,18 @@ int main(int argc, char *argv[])
     );
     argList::addBoolOption
     (
-        "auto-origin",
-        "Use bounding box centre as origin for rotations"
+        "auto-centre",
+        "Use bounding box centre as centre for rotations"
     );
     argList::addOption
     (
-        "origin",
+        "centre",
         "point",
-        "Use specified <point> as origin for rotations"
+        "Use specified <point> as centre for rotations"
     );
+    argList::addOptionCompat("auto-centre", {"auto-origin", 2206});
+    argList::addOptionCompat("centre", {"origin", 2206});
+
     argList::addOption
     (
         "rotate",
@@ -352,18 +355,18 @@ int main(int argc, char *argv[])
         points += v;
     }
 
-    vector origin;
-    bool useOrigin = args.readIfPresent("origin", origin);
-    if (args.found("auto-origin") && !useOrigin)
+    vector rotationCentre;
+    bool useRotationCentre = args.readIfPresent("centre", rotationCentre);
+    if (args.found("auto-centre") && !useRotationCentre)
     {
-        useOrigin = true;
-        origin = boundBox(points).centre();
+        useRotationCentre = true;
+        rotationCentre = boundBox(points).centre();
     }
 
-    if (useOrigin)
+    if (useRotationCentre)
     {
-        Info<< "Set origin for rotations to " << origin << endl;
-        points -= origin;
+        Info<< "Set centre of rotation to " << rotationCentre << endl;
+        points -= rotationCentre;
     }
 
 
@@ -455,15 +458,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);
 

bin/foamCleanPath

@@ -7,7 +7,7 @@
 #    \\/     M anipulation  |
 #------------------------------------------------------------------------------
 #     Copyright (C) 2011-2016 OpenFOAM Foundation
-#     Copyright (C) 2017-2021 OpenCFD Ltd.
+#     Copyright (C) 2017-2022 OpenCFD Ltd.
 #------------------------------------------------------------------------------
 # License
 #     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -48,26 +48,29 @@
 #     - Similarly for c-shell
 #       eval `foamCleanPath -csh-path dir1:dir2`
 #
-# 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.
+# 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.
 #
 #------------------------------------------------------------------------------
 printHelp() {
     cat<<USAGE
 
-Usage: foamCleanPath [OPTION] ENV [filter] .. [filter]
+Usage: foamCleanPath [OPTION] ENVNAME [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 (DYLD_LIBRARY_PATH on Darwin)
+                    and LD_LIBRARY_PATH (FOAM_LD_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
-                        or DYLD_LIBRARY_PATH on Darwin
+                        (FOAM_LD_LIBRARY_PATH on Darwin)
   -debug            Print debug information to stderr
   -strip            Remove inaccessible directories
   -verbose          Report some progress (input, output, ...)
@@ -81,7 +84,7 @@ Prints its argument (which should be a ':' separated list) cleansed from
 Exit status
     0  on success
     1  general error
-    2  initial value of 'path' is empty
+    2  initial value of ENVNAME is empty
 
 USAGE
     exit 0  # A clean exit
@@ -103,63 +106,69 @@ die()
 #-------------------------------------------------------------------------------
 
 # Input and outputs
-unset dirList shellOutput shellFlavour
+unset dirList shellFlavour shellOutput
 unset optDebug optEnvName optStrip optVerbose
 
 # Parse options
 while [ "$#" -gt 0 ]
 do
     case "$1" in
-    -h | -help*)
+    (-h | -help*)
         printHelp
         ;;
 
-    -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
+    (-csh-lib | -csh-path | -sh-lib | -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'
             ;;
-         esac
-         optEnvName="$name"     # Use for input evaluation
-         ;;
-
-    -csh-path | -sh-path)
-        shellFlavour="$1"
-        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
         ;;
 
-    -csh=* | -sh=* | -csh-env=* | -sh-env=*)
+    (-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=*)
         shellFlavour="$1"
         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
+        # 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
         shellOutput="$name"     # Use for output
@@ -167,16 +176,16 @@ do
         case "$1" in (*-env=*) optEnvName="$name";; esac
         ;;
 
-    -debug)
+    (-debug)
         optDebug=true
         ;;
-    -strip)
+    (-strip)
         optStrip=true
         ;;
-    -verbose)
+    (-verbose)
         optVerbose=true
         ;;
-    *)
+    (*)
         break
         ;;
     esac