CONFIG: bump patch level (240625)

- the old logic relied on the presence/absence of a comma to decide
  whether to parse as fixed or free format. This logic is faulty when
  handling the final (trailing) entry and will generally lead to the
  last field being truncated when read in.
  Now the caller decides on fixed vs free.

FIX: inconsistent Nastran surface output format

- use FREE format by default. Previously had an odd mix of SHORT
  format when created without options and LONG format (as default)
  when created with format options.

META-INFO/api-info

@@ -1,2 +1,2 @@
 api=2212
-patch=230110
+patch=240625

applications/solvers/compressible/rhoPimpleFoam/Make/options

@@ -9,7 +9,7 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/regionFaModels\lnInclude
+    -I$(LIB_SRC)/regionFaModels/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \

applications/solvers/compressible/rhoSimpleFoam/Make/options

@@ -7,7 +7,7 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/regionFaModels\lnInclude
+    -I$(LIB_SRC)/regionFaModels/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \

applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options

@@ -10,7 +10,7 @@ EXE_INC = \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-    -I$(LIB_SRC)/regionFaModels\lnInclude
+    -I$(LIB_SRC)/regionFaModels/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \

applications/solvers/heatTransfer/buoyantSimpleFoam/Make/options

@@ -7,7 +7,7 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/regionFaModels\lnInclude
+    -I$(LIB_SRC)/regionFaModels/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \

applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options

@@ -19,7 +19,7 @@ EXE_INC = \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
-    -I$(LIB_SRC)/regionFaModels\lnInclude
+    -I$(LIB_SRC)/regionFaModels/lnInclude
 
 
 EXE_LIBS = \

applications/solvers/incompressible/pimpleFoam/Make/options

@@ -8,7 +8,7 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-    -I$(LIB_SRC)/regionFaModels\lnInclude
+    -I$(LIB_SRC)/regionFaModels/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/CGALTriangulation3DKernel.H

@@ -30,8 +30,8 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef CGALTriangulation3DKernel_H
-#define CGALTriangulation3DKernel_H
+#ifndef Foam_CGALTriangulation3DKernel_H
+#define Foam_CGALTriangulation3DKernel_H
 
 // Silence boost bind deprecation warnings (before CGAL-5.2.1)
 #include "CGAL/version.h"
@@ -54,9 +54,19 @@ Description
 //    #include "CGAL/Robust_circumcenter_traits_3.h"
 //    typedef CGAL::Robust_circumcenter_traits_3<baseK>           K;
 
+#if defined(CGAL_VERSION_NR) && (CGAL_VERSION_NR < 1050500000)
+    // Prior to CGAL-5.5
+
     #include "CGAL/Robust_circumcenter_filtered_traits_3.h"
     typedef CGAL::Robust_circumcenter_filtered_traits_3<baseK>  K;
 
+#else
+
+    #include "CGAL/Robust_weighted_circumcenter_filtered_traits_3.h"
+    typedef CGAL::Robust_weighted_circumcenter_filtered_traits_3<baseK> K;
+
+#endif
+
 #else
 
     // Very robust but expensive kernel

applications/utilities/preProcessing/viewFactorsGen/Make/options

@@ -4,7 +4,6 @@ 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 \
@@ -12,7 +11,6 @@ EXE_INC = \
 
 
 EXE_LIBS = \
-    /* ${CGAL_LIBS} */ \
     -lfiniteVolume \
     -lsurfMesh \
     -lmeshTools \

etc/caseDicts/postProcessing/solvers/scalarTransport/scalarTransport

@@ -18,6 +18,8 @@ Description
 type            scalarTransport;
 libs            ("libsolverFunctionObjects.so");
 
+writeControl    writeTime;
+
 field           s;
 schemesField    s;
 D               1e-09;

src/OpenFOAM/db/dictionary/functionEntries/evalEntry/evalEntry.C

@@ -174,7 +174,7 @@ Foam::tokenList Foam::functionEntries::evalEntry::evaluate
         result.writeField(toks);
     }
 
-    return std::move(toks);
+    return tokenList(std::move(toks.tokens()));
 }
 
 

src/OpenFOAM/expressions/exprDriver/exprDriverFields.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -41,9 +41,11 @@ void Foam::expressions::exprDriver::fill_random
 {
     if (seed <= 0)
     {
-        if (timeStatePtr_)
+        const TimeState* ts = this->timeState();
+
+        if (ts)
         {
-            seed = (timeStatePtr_->timeIndex() - seed);
+            seed = (ts->timeIndex() - seed);
         }
         else
         {

src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.C

@@ -152,10 +152,11 @@ slicedBoundaryField
                 new SlicedPatchField<Type>
                 (
                     mesh.boundary()[patchi],
-                    DimensionedField<Type, GeoMesh>::null(),
-                    bField[patchi]
+                    DimensionedField<Type, GeoMesh>::null()
                 )
             );
+
+            bf[patchi].UList<Type>::shallowCopy(bField[patchi]);
         }
     }
 

src/OpenFOAM/matrices/Matrix/Matrix.C

@@ -201,7 +201,8 @@ inline Foam::Matrix<Form, Type>::Matrix
 )
 :
     mRows_(Mb.m()),
-    nCols_(Mb.n())
+    nCols_(Mb.n()),
+    v_(nullptr)
 {
     doAlloc();
 
@@ -223,7 +224,8 @@ inline Foam::Matrix<Form, Type>::Matrix
 )
 :
     mRows_(Mb.m()),
-    nCols_(Mb.n())
+    nCols_(Mb.n()),
+    v_(nullptr)
 {
     doAlloc();
 

src/OpenFOAM/parallel/globalIndex/globalIndex.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -489,12 +489,12 @@ public:
 
         //- Collect indirect data in processor order on master
         //  Handles contiguous/non-contiguous data, skips empty fields.
-        template<class Type, class Addr>
+        template<class ProcIDsContainer, class Type, class Addr>
         static void gather
         (
             const labelUList& offsets,  //!< offsets (master only)
             const label comm,           //!< communicator
-            const UList<int>& procIDs,
+            const ProcIDsContainer& procIDs,
             const IndirectListBase<Type, Addr>& fld,
             List<Type>& allFld,         //! output field (master only)
             const int tag = UPstream::msgType(),

src/OpenFOAM/parallel/globalIndex/globalIndexTemplates.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2017 OpenFOAM Foundation
-    Copyright (C) 2019-2022 OpenCFD Ltd.
+    Copyright (C) 2019-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -248,12 +248,12 @@ void Foam::globalIndex::gather
 }
 
 
-template<class Type, class Addr>
+template<class ProcIDsContainer, class Type, class Addr>
 void Foam::globalIndex::gather
 (
     const labelUList& off,  // needed on master only
     const label comm,
-    const UList<int>& procIDs,
+    const ProcIDsContainer& procIDs,
     const IndirectListBase<Type, Addr>& fld,
     List<Type>& allFld,
     const int tag,
@@ -368,7 +368,7 @@ void Foam::globalIndex::gather
         (
             offsets_,  // needed on master only
             comm,
-            UPstream::procID(comm),
+            UPstream::allProcs(comm),   // All communicator ranks
             sendData,
             allData,
             tag,
@@ -404,7 +404,7 @@ void Foam::globalIndex::gather
         (
             offsets_,  // needed on master only
             comm,
-            UPstream::procID(comm),
+            UPstream::allProcs(comm),   // All communicator ranks
             sendData,
             allData,
             tag,
@@ -622,7 +622,7 @@ void Foam::globalIndex::mpiGather
             (
                 offsets_,  // needed on master only
                 comm,
-                UPstream::procID(comm),
+                UPstream::allProcs(comm),   // All communicator ranks
                 sendData,
                 allData,
                 tag,
@@ -967,7 +967,7 @@ void Foam::globalIndex::scatter
         (
             offsets_,  // needed on master only
             comm,
-            UPstream::procID(comm),
+            UPstream::allProcs(comm),  // All communicator ranks
             allData,
             localData,
             tag,

src/OpenFOAM/primitives/Tensor/TensorI.H

@@ -409,6 +409,10 @@ inline Foam::Tensor<Cmpt> Foam::Tensor<Cmpt>::T() const
 
 
 template<class Cmpt>
+#if defined(__GNUC__) && !defined(__clang__)
+// Workaround for gcc (11+) that fails to handle tensor dot vector
+__attribute__((optimize("no-tree-vectorize")))
+#endif
 inline Foam::Tensor<Cmpt>
 Foam::Tensor<Cmpt>::inner(const Tensor<Cmpt>& t2) const
 {
@@ -432,6 +436,10 @@ Foam::Tensor<Cmpt>::inner(const Tensor<Cmpt>& t2) const
 
 
 template<class Cmpt>
+#if defined(__GNUC__) && !defined(__clang__)
+// Workaround for gcc (11+) that fails to handle tensor dot vector
+__attribute__((optimize("no-tree-vectorize")))
+#endif
 inline Foam::Tensor<Cmpt>
 Foam::Tensor<Cmpt>::schur(const Tensor<Cmpt>& t2) const
 {
@@ -867,6 +875,10 @@ operator&(const Tensor<Cmpt>& t1, const Tensor<Cmpt>& t2)
 
 //- Inner-product of a SphericalTensor and a Tensor
 template<class Cmpt>
+#if defined(__GNUC__) && !defined(__clang__)
+// Workaround for gcc (11+) that fails to handle tensor dot vector
+__attribute__((optimize("no-tree-vectorize")))
+#endif
 inline Tensor<Cmpt>
 operator&(const SphericalTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)
 {
@@ -881,6 +893,10 @@ operator&(const SphericalTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)
 
 //- Inner-product of a Tensor and a SphericalTensor
 template<class Cmpt>
+#if defined(__GNUC__) && !defined(__clang__)
+// Workaround for gcc (11+) that fails to handle tensor dot vector
+__attribute__((optimize("no-tree-vectorize")))
+#endif
 inline Tensor<Cmpt>
 operator&(const Tensor<Cmpt>& t1, const SphericalTensor<Cmpt>& st2)
 {
@@ -895,6 +911,10 @@ operator&(const Tensor<Cmpt>& t1, const SphericalTensor<Cmpt>& st2)
 
 //- Inner-product of a SymmTensor and a Tensor
 template<class Cmpt>
+#if defined(__GNUC__) && !defined(__clang__)
+// Workaround for gcc (11+) that fails to handle tensor dot vector
+__attribute__((optimize("no-tree-vectorize")))
+#endif
 inline Tensor<Cmpt>
 operator&(const SymmTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)
 {
@@ -917,6 +937,10 @@ operator&(const SymmTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)
 
 //- Inner-product of a Tensor and a SymmTensor
 template<class Cmpt>
+#if defined(__GNUC__) && !defined(__clang__)
+// Workaround for gcc (11+) that fails to handle tensor dot vector
+__attribute__((optimize("no-tree-vectorize")))
+#endif
 inline Tensor<Cmpt>
 operator&(const Tensor<Cmpt>& t1, const SymmTensor<Cmpt>& st2)
 {
@@ -957,6 +981,10 @@ operator&(const Tensor<Cmpt>& t, const Vector<Cmpt>& v)
 
 //- Inner-product of a Vector and a Tensor
 template<class Cmpt>
+#if defined(__GNUC__) && !defined(__clang__)
+// Workaround for gcc (11+) that fails to handle tensor dot vector
+__attribute__((optimize("no-tree-vectorize")))
+#endif
 inline Vector<Cmpt>
 operator&(const Vector<Cmpt>& v, const Tensor<Cmpt>& t)
 {

src/OpenFOAM/primitives/strings/word/word.C

@@ -29,6 +29,7 @@ License
 #include "word.H"
 #include "debug.H"
 #include <cctype>
+#include <cstdint>
 #include <sstream>
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

src/Pstream/mpi/UPstream.C

@@ -641,7 +641,16 @@ void Foam::UPstream::freePstreamCommunicator(const label communicator)
     }
 
     // Not touching the first two communicators (SELF, WORLD)
-    if (communicator > 1)
+    // or anything out-of bounds.
+    //
+    // No UPstream communicator indices when MPI is initialized outside
+    // of OpenFOAM - thus needs a bounds check too!
+
+    if
+    (
+        communicator > 1
+     && (communicator < PstreamGlobals::MPICommunicators_.size())
+    )
     {
         if (MPI_COMM_NULL != PstreamGlobals::MPICommunicators_[communicator])
         {

src/TurbulenceModels/compressible/turbulentFluidThermoModels/turbulentFluidThermoModels.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
-    Copyright (C) 2022 OpenCFD Ltd.
+    Copyright (C) 2022-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -131,6 +131,9 @@ makeLESModel(dynamicKEqn);
 #include "dynamicLagrangian.H"
 makeLESModel(dynamicLagrangian);
 
+#include "sigma.H"
+makeLESModel(sigma);
+
 #include "SpalartAllmarasDES.H"
 makeLESModel(SpalartAllmarasDES);
 

src/TurbulenceModels/turbulenceModels/laminar/generalizedNewtonian/generalizedNewtonianViscosityModels/Casson/Casson.C

@@ -79,8 +79,8 @@ bool Foam::laminarModels::generalizedNewtonianViscosityModels::Casson::read
 
     coeffs.readEntry("m", m_);
     coeffs.readEntry("tau0", tau0_);
-    coeffs.readEntry("nuMin_", nuMin_);
-    coeffs.readEntry("nuMax_", nuMax_);
+    coeffs.readEntry("nuMin", nuMin_);
+    coeffs.readEntry("nuMax", nuMax_);
 
     return true;
 }

src/fileFormats/nastran/NASCore.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2022 OpenCFD Ltd.
+    Copyright (C) 2017-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -122,7 +122,8 @@ std::string Foam::fileFormats::NASCore::nextNasField
 (
     const std::string& str,
     std::string::size_type& pos,
-    std::string::size_type len
+    const std::string::size_type width,
+    const bool free_format
 )
 {
     const auto beg = pos;
@@ -130,15 +131,23 @@ std::string Foam::fileFormats::NASCore::nextNasField
 
     if (end == std::string::npos)
     {
-        pos = beg + len;    // Continue after field width
+        if (free_format)
+        {
+            // Nothing left
+            pos = str.size();
+            return str.substr(beg);
+        }
+
+        // Fixed format - continue after field width
+        pos = beg + width;
+        return str.substr(beg, width);
     }
     else
     {
-        len = (end - beg);  // Efffective width
-        pos = end + 1;      // Continue after comma
+        // Free format - continue after comma
+        pos = end + 1;
+        return str.substr(beg, (end - beg));
     }
-
-    return str.substr(beg, len);
 }
 
 
@@ -235,8 +244,8 @@ void Foam::fileFormats::NASCore::writeCoord
     // 2 ID   : point ID - requires starting index of 1
     // 3 CP   : coordinate system ID                (blank)
     // 4 X1   : point x coordinate
-    // 5 X2   : point x coordinate
-    // 6 X3   : point x coordinate
+    // 5 X2   : point y coordinate
+    // 6 X3   : point z coordinate
     // 7 CD   : coordinate system for displacements (blank)
     // 8 PS   : single point constraints            (blank)
     // 9 SEID : super-element ID

src/fileFormats/nastran/NASCore.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2022 OpenCFD Ltd.
+    Copyright (C) 2017-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -48,7 +48,6 @@ SourceFiles
 
 namespace Foam
 {
-
 namespace fileFormats
 {
 
@@ -74,18 +73,18 @@ public:
         //- Output load format
         enum loadFormat
         {
-            PLOAD2,
-            PLOAD4
+            PLOAD2,     //!< Face load (eg, pressure)
+            PLOAD4      //!< Vertex load
         };
 
-        //- Selection names for the NASTRAN file field formats
+        //- Selection names for the NASTRAN load formats
         static const Enum<loadFormat> loadFormatNames;
 
 
     // Constructors
 
         //- Default construct
-        NASCore() = default;
+        NASCore() noexcept = default;
 
 
     // Public Static Member Functions
@@ -93,18 +92,20 @@ public:
         //- Extract numbers from things like "-2.358-8" (same as "-2.358e-8")
         static scalar readNasScalar(const std::string& str);
 
-        //- A string::substr() to handle fixed-format and free-format NASTRAN.
-        //  Returns the substr to the next comma (if found) or the given length
-        //
-        //  \param str The string to extract from
-        //  \param pos On input, the position of the first character of the
-        //      substring. On output, advances to the next position to use.
-        //  \param len The fixed-format length to use if a comma is not found.
+        //- A std::string::substr() variant to handle fixed-format and
+        //- free-format NASTRAN.
+        //  Returns the substr until the next comma (if found)
+        //  or the given fixed width
         static std::string nextNasField
         (
+            //! The string to extract from
             const std::string& str,
+            //! [in,out] The parse position within \p str
             std::string::size_type& pos,
-            std::string::size_type len
+            //! The fixed-format width to use (if comma is not found)
+            const std::string::size_type width,
+            //! The input is known to be free-format
+            const bool free_format = false
         );
 
 

src/fileFormats/stl/STLCore.C

@@ -135,8 +135,9 @@ int Foam::fileFormats::STLCore::detectBinaryHeader
 
         bad =
             (
-                nTris < int(dataFileSize - STLHeaderSize)/50
-             || nTris > int(dataFileSize - STLHeaderSize)/25
+                dataFileSize < STLHeaderSize
+             || nTris < (dataFileSize - STLHeaderSize)/50
+             || nTris > (dataFileSize - STLHeaderSize)/25
             );
     }
 
@@ -208,8 +209,9 @@ Foam::fileFormats::STLCore::readBinaryHeader
 
         bad =
             (
-                nTris < int(dataFileSize - STLHeaderSize)/50
-             || nTris > int(dataFileSize - STLHeaderSize)/25
+                dataFileSize < STLHeaderSize
+             || nTris < (dataFileSize - STLHeaderSize)/50
+             || nTris > (dataFileSize - STLHeaderSize)/25
             );
     }
 

src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -180,12 +180,26 @@ Foam::faBoundaryMesh::faBoundaryMesh
 
 void Foam::faBoundaryMesh::calcGeometry()
 {
-    // processorFaPatch geometry triggers calculation of pointNormals.
+    // processor initGeometry send/recv the following:
+    //   - edgeCentres() : faMesh::edgeCentres()
+    //   - edgeLengths() : faMesh::Le()
+    //   - edgeFaceCentres() : faMesh::areaCentres()
+    //
+    // faMesh::Le() has its own point-to-point communication (OK) but
+    // triggers either/or edgeAreaNormals(), pointAreaNormals()
+    // with their own communication that can block.
+
     // This uses parallel comms and hence will not be trigggered
     // on processors that do not have a processorFaPatch so instead
     // force construction.
+
+    (void)mesh_.edgeAreaNormals();
     (void)mesh_.pointAreaNormals();
 
+    (void)mesh_.areaCentres();
+    (void)mesh_.faceAreaNormals();
+
+
     PstreamBuffers pBufs(Pstream::defaultCommsType);
 
     if
@@ -773,12 +787,15 @@ bool Foam::faBoundaryMesh::checkDefinition(const bool report) const
 
 void Foam::faBoundaryMesh::movePoints(const pointField& p)
 {
-    // processorFaPatch geometry triggers calculation of pointNormals.
-    // This uses parallel comms and hence will not be trigggered
-    // on processors that do not have a processorFaPatch so instead
-    // force construction.
+    // See comments in calcGeometry()
+
+    (void)mesh_.edgeAreaNormals();
     (void)mesh_.pointAreaNormals();
 
+    (void)mesh_.areaCentres();
+    (void)mesh_.faceAreaNormals();
+
+
     PstreamBuffers pBufs(Pstream::defaultCommsType);
 
     if

src/finiteArea/faMesh/faMesh.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2020-2022 OpenCFD Ltd.
+    Copyright (C) 2020-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -277,6 +277,27 @@ void Foam::faMesh::clearOut() const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+void Foam::faMesh::syncGeom()
+{
+    if (UPstream::parRun())
+    {
+        // areaCentres()
+        if (faceCentresPtr_)
+        {
+            faceCentresPtr_->boundaryFieldRef()
+                .evaluateCoupled<processorFaPatch>();
+        }
+
+        // faceAreaNormals()
+        if (faceAreaNormalsPtr_)
+        {
+            faceAreaNormalsPtr_->boundaryFieldRef()
+                .evaluateCoupled<processorFaPatch>();
+        }
+    }
+}
+
+
 bool Foam::faMesh::init(const bool doInit)
 {
     if (doInit)
@@ -296,18 +317,7 @@ bool Foam::faMesh::init(const bool doInit)
     // Calculate the geometry for the patches (transformation tensors etc.)
     boundary_.calcGeometry();
 
-    // Ensure processor/processor information is properly synchronised
-    if (Pstream::parRun())
-    {
-        const_cast<areaVectorField&>(areaCentres()).boundaryFieldRef()
-            .evaluateCoupled<processorFaPatch>();
-
-        // This roughly corresponds to what OpenFOAM-v2112 (and earlier) had,
-        // but should nominally be unnecessary.
-        //
-        /// const_cast<areaVectorField&>(faceAreaNormals()).boundaryFieldRef()
-        ///     .evaluateCoupled<processorFaPatch>();
-    }
+    syncGeom();
 
     return false;
 }
@@ -989,7 +999,6 @@ bool Foam::faMesh::movePoints()
 
     clearGeomNotAreas();
 
-    // To satisfy the motion interface for MeshObject, const cast is needed
     if (patchPtr_)
     {
         patchPtr_->movePoints(newPoints);
@@ -1003,6 +1012,8 @@ bool Foam::faMesh::movePoints()
 
     // Note: Fluxes were dummy?
 
+    syncGeom();
+
     return true;
 }
 

src/finiteArea/faMesh/faMesh.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2021-2022 OpenCFD Ltd.
+    Copyright (C) 2021-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -617,6 +617,10 @@ public:
         //- Initialise non-demand-driven data etc
         bool init(const bool doInit);
 
+        //- Processor/processor synchronisation for geometry fields.
+        //  Largely internal use only (slightly hacky).
+        void syncGeom();
+
 
     // Database
 

src/finiteArea/faMesh/faMeshDemandDrivenData.C

@@ -898,6 +898,12 @@ void Foam::faMesh::calcFaceCentres() const
             }
         }
     }
+
+    // Parallel consistency, exchange on processor patches
+    if (UPstream::parRun())
+    {
+        centres.boundaryFieldRef().evaluateCoupled<processorFaPatch>();
+    }
 }
 
 
@@ -1110,6 +1116,12 @@ void Foam::faMesh::calcFaceAreaNormals() const
                 = edgeNormalsBoundary[patchi];
         }
     }
+
+    // Parallel consistency, exchange on processor patches
+    if (UPstream::parRun())
+    {
+        faceNormals.boundaryFieldRef().evaluateCoupled<processorFaPatch>();
+    }
 }
 
 

src/finiteArea/faMesh/faMeshTools/faMeshTools.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2022 OpenCFD Ltd.
+    Copyright (C) 2015-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -61,6 +61,8 @@ void Foam::faMeshTools::forceDemandDriven(faMesh& mesh)
     (void)mesh.pointAreaNormals();
     (void)mesh.faceCurvatures();
     (void)mesh.edgeTransformTensors();
+
+    mesh.syncGeom();
 }
 
 

src/functionObjects/field/expressions/fvExpressionField.C

@@ -389,7 +389,7 @@ bool Foam::functionObjects::fvExpressionField::read(const dictionary& dict)
     }
 
     autowrite_ = dict.getOrDefault("autowrite", false);
-    store_ = dict.getOrDefault("autowrite", true);
+    store_ = dict.getOrDefault("store", true);
 
     // "dimensions" is optional
     dimensions_.clear();

src/functionObjects/utilities/ensightWrite/ensightWriteUpdate.C

@@ -103,7 +103,7 @@ bool Foam::functionObjects::ensightWrite::readSelection(const dictionary& dict)
     selectFields_.uniq();
 
     blockFields_.clear();
-    dict.readIfPresent("blockFields", blockFields_);
+    dict.readIfPresent("excludeFields", blockFields_);
     blockFields_.uniq();
 
     // Actions to define selection

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectedParticleDistributionInjection/InjectedParticleDistributionInjection.C

@@ -179,7 +179,7 @@ void Foam::InjectedParticleDistributionInjection<CloudType>::initialise()
                 sumPow3 += pow3(diameters[particlei]);
             }
 
-            const scalar volume = sumPow3*mathematical::pi/16.0;
+            const scalar volume = sumPow3*mathematical::pi/6.0;
             sumVolume += volume;
             volumeFlowRate_[injectori] = volume/dTime;
 

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectedParticleInjection/InjectedParticleInjection.C

@@ -119,7 +119,7 @@ void Foam::InjectedParticleInjection<CloudType>::initialise()
     scalar sumVolume = 0;
     forAll(volume, i)
     {
-        scalar vol = pow3(diameter_[i])*mathematical::pi/16.0;
+        scalar vol = pow3(diameter_[i])*mathematical::pi/6.0;
         volume[i] = vol;
         sumVolume += vol;
     }

src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2019-2022 OpenCFD Ltd.
+    Copyright (C) 2019-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -62,6 +62,53 @@ Foam::ThermoSurfaceFilm<CloudType>::ThermoSurfaceFilm
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class CloudType>
+template<class filmType>
+void Foam::ThermoSurfaceFilm<CloudType>::absorbInteraction
+(
+    filmType& film,
+    const parcelType& p,
+    const polyPatch& pp,
+    const label facei,
+    const scalar mass,
+    bool& keepParticle
+)
+{
+    DebugInfo<< "Parcel " << p.origId() << " absorbInteraction" << endl;
+
+    // Patch face normal
+    const vector& nf = pp.faceNormals()[facei];
+
+    // Patch velocity
+    const vector& Up = this->owner().U().boundaryField()[pp.index()][facei];
+
+    // Relative parcel velocity
+    const vector Urel(p.U() - Up);
+
+    // Parcel normal velocity
+    const vector Un(nf*(Urel & nf));
+
+    // Parcel tangential velocity
+    const vector Ut(Urel - Un);
+
+    film.addSources
+    (
+        pp.index(),
+        facei,
+        mass,                             // mass
+        mass*Ut,                          // tangential momentum
+        mass*mag(Un),                     // impingement pressure
+        mass*p.hs()                       // energy
+    );
+
+    this->nParcelsTransferred()++;
+
+    this->totalMassTransferred() += mass;
+
+    keepParticle = false;
+}
+
+
 template<class CloudType>
 bool Foam::ThermoSurfaceFilm<CloudType>::transferParcel
 (

src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2021-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -143,6 +143,21 @@ public:
 
     // Member Functions
 
+        // Interaction models
+
+            //- Absorb parcel into film
+            template<class filmType>
+            void absorbInteraction
+            (
+                filmType&,
+                const parcelType& p,
+                const polyPatch& pp,
+                const label facei,
+                const scalar mass,
+                bool& keepParticle
+            );
+
+
         // Evaluation
 
             //- Transfer parcel from cloud to surface film

src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C

@@ -1368,6 +1368,7 @@ Foam::label Foam::snappyRefineDriver::refinementInterfaceRefine
                                 (
                                     face2i != facei
                                  && surfaceIndex[face2i] != -1
+                                 && cutter.faceLevel(face2i) > cLevel
                                 )
                                 {
                                     // Get outwards pointing normal

src/meshTools/AMIInterpolation/AMIInterpolation/faceAreaWeightAMI/faceAreaWeightAMI.C

@@ -269,7 +269,8 @@ bool Foam::faceAreaWeightAMI::setNextFaces
         return false;
     }
 
-    const labelList& srcNbrFaces = this->srcPatch().faceFaces()[srcFacei];
+    const auto& srcPatch = this->srcPatch();
+    const labelList& srcNbrFaces = srcPatch.faceFaces()[srcFacei];
 
     // Initialise tgtFacei
     tgtFacei = -1;
@@ -360,6 +361,7 @@ bool Foam::faceAreaWeightAMI::setNextFaces
     {
         FatalErrorInFunction
             << "Unable to set target face for source face " << srcFacei
+            << " with centre: " << srcPatch.faceCentres()[srcFacei]
             << abort(FatalError);
     }
 

src/meshTools/edgeMesh/edgeFormats/nastran/NASedgeFormat.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2017 OpenCFD Ltd.
+    Copyright (C) 2017-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -62,13 +62,17 @@ bool Foam::fileFormats::NASedgeFormat::read
 
     while (is.good())
     {
-        string::size_type linei = 0;  // parsing position within current line
         string line;
         is.getLine(line);
 
-        if (line.empty() || line[0] == '$')
+        if (line.empty())
         {
-            continue; // Skip empty or comment
+            continue;  // Ignore empty
+        }
+        else if (line[0] == '$')
+        {
+            // Ignore comment
+            continue;
         }
 
         // Check if character 72 is continuation
@@ -94,38 +98,66 @@ bool Foam::fileFormats::NASedgeFormat::read
         }
 
 
+        // Parsing position within current line
+        std::string::size_type linei = 0;
+
+        // Is free format if line contains a comma
+        const bool freeFormat = line.contains(',');
+
         // First word (column 0-8)
         const word cmd(word::validate(nextNasField(line, linei, 8)));
 
         if (cmd == "CBEAM" || cmd == "CROD")
         {
-            // discard elementId (8-16)
-            (void) nextNasField(line, linei, 8); // 8-16
-            // discard groupId (16-24)
-            (void) nextNasField(line, linei, 8); // 16-24
+            // Fixed format:
+            //  8-16 : element id
+            // 16-24 : group id
+            // 24-32 : vertex
+            // 32-40 : vertex
 
-            label a = readLabel(nextNasField(line, linei, 8)); // 24-32
-            label b = readLabel(nextNasField(line, linei, 8)); // 32-40
+            // discard elementId
+            (void) nextNasField(line, linei, 8, freeFormat);
+            // discard groupId
+            (void) nextNasField(line, linei, 8, freeFormat);
+
+            label a = readLabel(nextNasField(line, linei, 8, freeFormat));
+            label b = readLabel(nextNasField(line, linei, 8, freeFormat));
 
             dynEdges.append(edge(a,b));
         }
         else if (cmd == "PLOTEL")
         {
-            // discard elementId (8-16)
-            (void) nextNasField(line, linei, 8); // 8-16
+            // Fixed format:
+            //  8-16 : element id
+            // 16-24 : vertex
+            // 24-32 : vertex
+            // 32-40 : vertex
 
-            label a = readLabel(nextNasField(line, linei, 8)); // 16-24
-            label b = readLabel(nextNasField(line, linei, 8)); // 24-32
+            // discard elementId (8-16)
+            (void) nextNasField(line, linei, 8, freeFormat);
+
+            label a = readLabel(nextNasField(line, linei, 8, freeFormat));
+            label b = readLabel(nextNasField(line, linei, 8, freeFormat));
 
             dynEdges.append(edge(a,b));
         }
         else if (cmd == "GRID")
         {
-            label index = readLabel(nextNasField(line, linei, 8)); // 8-16
-            (void) nextNasField(line, linei, 8); // 16-24
-            scalar x = readNasScalar(nextNasField(line, linei, 8)); // 24-32
-            scalar y = readNasScalar(nextNasField(line, linei, 8)); // 32-40
-            scalar z = readNasScalar(nextNasField(line, linei, 8)); // 40-48
+            // Fixed (short) format:
+            //  8-16 : point id
+            // 16-24 : coordinate system (unsupported)
+            // 24-32 : point x coordinate
+            // 32-40 : point y coordinate
+            // 40-48 : point z coordinate
+            // 48-56 : displacement coordinate system (optional, unsupported)
+            // 56-64 : single point constraints (optional, unsupported)
+            // 64-70 : super-element id (optional, unsupported)
+
+            label index = readLabel(nextNasField(line, linei, 8, freeFormat));
+            (void) nextNasField(line, linei, 8, freeFormat);
+            scalar x = readNasScalar(nextNasField(line, linei, 8, freeFormat));
+            scalar y = readNasScalar(nextNasField(line, linei, 8, freeFormat));
+            scalar z = readNasScalar(nextNasField(line, linei, 8, freeFormat));
 
             pointId.append(index);
             dynPoints.append(point(x, y, z));
@@ -138,6 +170,8 @@ bool Foam::fileFormats::NASedgeFormat::read
             // GRID*      126   0 -5.55999875E+02 -5.68730474E+02
             // *         2.14897901E+02
 
+            // Cannot be long format and free format at the same time!
+
             label index = readLabel(nextNasField(line, linei, 16)); // 8-24
             (void) nextNasField(line, linei, 16); // 24-40
             scalar x = readNasScalar(nextNasField(line, linei, 16)); // 40-56

src/meshTools/output/vtk/topoSet/foamVtkWritePointSet.C

@@ -155,7 +155,7 @@ bool Foam::vtk::writePointSet
 
     if (parallel)
     {
-        vtk::writeListParallel(format(), mesh.points(), pointLabels);
+        vtk::writeListParallel(format.ref(), mesh.points(), pointLabels);
     }
     else
     {

src/optimisation/adjointOptimisation/adjoint/ATCModel/ATCUaGradU/ATCUaGradU.C

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2021 PCOpt/NTUA
-    Copyright (C) 2013-2021 FOSS GP
+    Copyright (C) 2007-2023 PCOpt/NTUA
+    Copyright (C) 2013-2023 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -86,10 +86,9 @@ void ATCUaGradU::addATC(fvVectorMatrix& UaEqn)
     if (extraConvection_ > 0)
     {
         // Implicit part added to increase diagonal dominance
-        // Note: Maybe this needs to be multiplied with the ATClimiter ??
-        UaEqn += extraConvection_*fvm::div(-phi, Ua);
+        UaEqn += ATClimiter_*extraConvection_*fvm::div(-phi, Ua);
 
-        // correct rhs due to implicitly augmenting the adjoint convection
+        // Correct rhs due to implicitly augmenting the adjoint convection
         ATC_ += extraConvection_*(fvc::grad(UaForATC(), "gradUaATC")().T() & U);
     }
 

src/optimisation/adjointOptimisation/adjoint/ATCModel/ATCstandard/ATCstandard.C

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2021 PCOpt/NTUA
-    Copyright (C) 2013-2021 FOSS GP
+    Copyright (C) 2007-2023 PCOpt/NTUA
+    Copyright (C) 2013-2023 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -90,10 +90,9 @@ void ATCstandard::addATC(fvVectorMatrix& UaEqn)
     if (extraConvection_ > 0)
     {
         // Implicit part added to increase diagonal dominance
-        // Note: Maybe this needs to be multiplied with the ATClimiter ??
-        UaEqn += extraConvection_*fvm::div(-phi, Ua);
+        UaEqn += ATClimiter_*extraConvection_*fvm::div(-phi, Ua);
 
-        // correct rhs due to implicitly augmenting the adjoint convection
+        // Correct rhs due to implicitly augmenting the adjoint convection
         ATC_ += extraConvection_*(fvc::grad(Ua, "gradUaATC")().T() & U);
     }
 

src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C

@@ -1949,7 +1949,7 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
 //        //globalBorderTris.gather
 //        //(
 //        //    UPstream::worldComm,
-//        //    UPstream::procID(Pstream::worldComm),
+//        //    UPstream::allProcs(UPstream::worldComm),
 //        //    globalBorderCentres
 //        //);
 //        pointField globalBorderCentres(allCentres);
@@ -1996,7 +1996,7 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
 //        //globalBorderTris.scatter
 //        //(
 //        //    UPstream::worldComm,
-//        //    UPstream::procID(Pstream::worldComm),
+//        //    UPstream::allProcs(UPstream::worldComm),
 //        //    isMasterPoint
 //        //);
 //        //boolList isMasterBorder(s.size(), false);
@@ -2094,7 +2094,7 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
             globalTris().gather
             (
                 UPstream::worldComm,
-                UPstream::procID(Pstream::worldComm),
+                UPstream::allProcs(UPstream::worldComm),
                 allCentres
             );
         }
@@ -2144,7 +2144,7 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
             globalTris().scatter
             (
                 UPstream::worldComm,
-                UPstream::procID(Pstream::worldComm),
+                UPstream::allProcs(UPstream::worldComm),
                 allDistribution,
                 distribution
             );

src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/alphatPhaseChangeWallFunction/alphatPhaseChangeWallFunctionFvPatchScalarField.C

@@ -75,7 +75,7 @@ alphatPhaseChangeWallFunctionFvPatchScalarField
 
     if (dict.found("mDotL"))
     {
-        dmdt_ = scalarField("mDotL", dict, p.size());
+        mDotL_ = scalarField("mDotL", dict, p.size());
     }
 }
 

src/renumber/renumberMethods/Make/files

@@ -5,6 +5,5 @@ randomRenumber/randomRenumber.C
 springRenumber/springRenumber.C
 structuredRenumber/structuredRenumber.C
 structuredRenumber/OppositeFaceCellWaveBase.C
-hpathRenumber/hpathRenumber.C
 
 LIB = $(FOAM_LIBBIN)/librenumberMethods

src/renumber/renumberMethods/hpathRenumber/hpathRenumber.C

@@ -1,976 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2023 Alon Zameret, Noam Manaker Morag
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "hpathRenumber.H"
-#include "addToRunTimeSelectionTable.H"
-
-#include "CircularBuffer.H"
-#include <iomanip>
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTypeNameAndDebug(hpathRenumber, 0);
-
-    addToRunTimeSelectionTable
-    (
-        renumberMethod,
-        hpathRenumber,
-        dictionary
-    );
-}
-
-
-// * * * * * * * * * * * * * * * * Local Details * * * * * * * * * * * * * * //
-
-
-/*---------------------------------------------------------------------------*\
-                 Class hpathRenumber::hpathFinder Declaration
-\*---------------------------------------------------------------------------*/
-
-namespace Foam
-{
-
-// Class hpathFinder Declaration
-class hpathRenumber::hpathFinder
-{
-    // Private Data
-
-        // The input mesh
-        const Foam::polyMesh& mesh;
-
-        // Counter for the number of renumbered cells
-        label nFoundCellCount;
-
-        // Marks cells that have been added to the renumbering as 'true'
-        Foam::bitSet bIsRenumbered;
-
-        // For every data structure, I explain what it is used for and which
-        // method is used to compute it
-
-        // For every cell, a list of all its point-neighbouring cells - getMeshGraph()
-        Foam::labelListList nMeshGraph;
-
-        // For each cell its 'layer index': - getLayers()
-        //      - cells in the same layer will have the same layer index
-        Foam::labelList nCellLayerIndex;
-
-        // For each cell its 'connected component index': - getConnectedComponents()
-        //      - cells in the same connected component will have the same connected component index
-        Foam::labelList nCellConnnectedComponentIndex;
-
-        // For each cell its face-neighbours in the connected component - getConnnectedComponentGraph()
-        Foam::labelListList nConnnectedComponentGraph;
-
-        // Marks cells that have already been by BFS - getStartingCellInConnnectedComponent()
-        Foam::bitSet bBFSFoundCell;
-
-        // For each cell its point distance from the connected components start cell   - reorderDistFromStart()
-        Foam::labelList nCellPointDistFromStart;
-        Foam::labelList nCellFaceDistFromStart;
-
-        // For each cell its DFS depth within the connected component     - findPath()
-        Foam::labelList nDFSCellDepth;
-
-        // For each cell its DFS parent within the connected component    - findPath()
-        Foam::labelList nDFSParentCell;
-
-
-    public: // Public methods
-
-        //- Constructor
-        explicit hpathFinder(const Foam::polyMesh& mesh);
-
-        // Member Functions
-
-        // Get Renumbering for the mesh
-        void getRenumbering
-        (
-            Foam::labelList& cellOrder,
-            bool bApplyLayerSeparation
-        );
-
-        // Returns the accuracy of the renumbering:
-        // Accuracy is defined as the percentage of consecutive cells that are also face-neighbours in the mesh
-        //  - Cells are face-neighbours if they have a common face
-        float getAccuracy(const Foam::labelList& cellOrder) const;
-
-
-    private: // Private methods
-
-        // Initialize data structures for later use
-        void initialize();
-
-        // Given a cell and a face index, find its neighbour through the face
-        // - If facing the boundary, returns -1
-        // - Otherwise, returns FaceOwner/FaceNeighbour[nFaceIdx], the one that's different from nCellIdx
-        label getNei(label nCellIdx, label nFaceIdx) const;
-
-        // Creates the 'Mesh-Graph': for every cell, a list of cells that are point-neighbours with it in the mesh
-        //  - Cells are point-neighbours if they have a common point
-        void getMeshGraph();
-
-        // Separates the mesh into layers: each cell has its layer saved in nCellLayerIndex
-        // Also returns for every layer a list of all cells in it
-        void getLayerSeparation(Foam::DynamicList<Foam::DynamicList<label>>& nCellsByLayer);
-
-        // For every connected component, find the deepest cell and choose it as a starting cell
-        // Returns a list with one starting cell per connected component
-        void getStartingCells(Foam::DynamicList<label>& nStartingCells) const;
-
-        // Once the starting cells have been found, this method does the actual layer separation
-        void getLayers(const Foam::labelList& nStartingCellList, Foam::DynamicList<Foam::DynamicList<label>>& nCellsByLayer);
-
-        // Renumber all cells in a layer
-        void solveLayer(Foam::labelList nCellsInLayer, Foam::labelList& cellOrder);
-
-        // Seperates cells into face-connected components
-        // This is done using a general DFS algorithm
-        void getConnectedComponents(const Foam::labelList& nCellList, Foam::DynamicList<Foam::DynamicList<label>>& nCellsByConnnectedComponent);
-
-        // Find an approximate H-path through a connected component
-        void solveConnectedComponent(const Foam::labelList& nCellsInConnectedComponent, Foam::labelList& cellOrder);
-
-        // Finds for each cell in the connected component a list of its face-neighbours within the component
-        void getConnnectedComponentGraph(const Foam::labelList& nCellsInConnectedComponent);
-
-        // Finds a starting cell within the connected component
-        label getStartingCellInConnnectedComponent(const Foam::labelList& nCellsInConnectedComponent);
-
-        // This method reorders the cells in the given connected component based on their distance from nStartCell
-        void reorderDistFromStart(label nStartCell, const Foam::labelList& nCellsInConnectedComponent);
-
-        // Finds an H-path within the connected component and returns it in nResultHpath
-        // H-path is guaranteed to start at nStartCell
-        void findPath(label nStartCell, Foam::labelList& cellOrder);
-
-        // Resets data structures for the cells that weren't found
-        void resetCells(Foam::labelList& nCellList);
-};
-
-} // End namespace Foam
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::hpathRenumber::hpathRenumber(const dictionary& dict)
-:
-    renumberMethod(dict),
-    applyLayerSeparation_
-    (
-        dict.optionalSubDict(typeName + "Coeffs")
-            .getOrDefault("layered", true)
-    )
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-Foam::labelList Foam::hpathRenumber::renumber
-(
-    const polyMesh& mesh,
-    const pointField& points
-) const
-{
-    Info<< nl
-        << "Starting Cell Renumbering: "
-        << mesh.nCells() << " Cells, "
-        << mesh.nFaces() << " Faces, "
-        << mesh.nPoints() << " Points" << nl << nl;
-
-    hpathRenumber::hpathFinder s(mesh);
-    Foam::labelList cellOrder;
-    s.getRenumbering(cellOrder, applyLayerSeparation_);
-
-    std::cout << std::endl << "Found path with accuracy of "
-        << std::fixed << std::setprecision(3)
-        << s.getAccuracy(cellOrder) << "%" << std::endl;
-
-    fprintf(stdout,"\n******************************************************\n\n\n");
-    return cellOrder;
-}
-
-
-// * * * * * * * * * * * * hpathRenumber::hpathFinder  * * * * * * * * * * * //
-
-// Public methods
-Foam::hpathRenumber::hpathFinder::hpathFinder(const polyMesh& mesh)
-:
-    mesh(mesh)
-{}
-
-
-void Foam::hpathRenumber::hpathFinder::getRenumbering
-(
-    Foam::labelList& cellOrder,
-    bool bApplyLayerSeparation
-)
-{
-    // Find a renumbering for the entire mesh
-
-    cellOrder.resize(mesh.nCells());
-    // Counter for how many cells we have added to the renumbering so far
-    nFoundCellCount = 0;
-
-    // Initialize the data structures:
-    Info<< "Initializing Data Structures" << endl;
-    initialize();
-
-    // Compute a graph to represent the mesh:
-    //      - Cells in the mesh will be connected in the graph if they have a *common point*
-    getMeshGraph();
-
-    Foam::DynamicList<Foam::DynamicList<label>> nCellsByLayer;
-    if (bApplyLayerSeparation)
-    {
-        getLayerSeparation(nCellsByLayer);
-    }
-    else
-    {
-        // If there is no layer separation, set the entire mesh as one 'layer'
-        Foam::DynamicList<label> nAllCells(Foam::identity(mesh.nCells()));
-        nCellsByLayer.append(nAllCells);
-    }
-
-    std::cout << "Beginning Hpath Computation" << std::endl;
-    // Find H-path for each layer separately
-    for (const Foam::labelList& nCellsInLayer : nCellsByLayer)
-    {
-        // solveLayer() will find a renumbering for all the cells in the layer
-        // Path will be appended into cellOrder
-        solveLayer(nCellsInLayer, cellOrder);
-    }
-}
-
-
-float Foam::hpathRenumber::hpathFinder::getAccuracy
-(
-    const Foam::labelList& cellOrder
-) const
-{
-    // Finding the number of "hits"
-    // - A hit is a pair of consecutive cells in the renumbering that are also face-neighbours in the mesh
-    // - Cells are face-neighbours if they have a common face
-    label nHitCnt = 0;
-    for (label nPathIdx = 0; nPathIdx < label(cellOrder.size()) - 1; nPathIdx++)
-    {
-        label nCurrCellIdx = cellOrder[nPathIdx];
-        label nNextCellIdx = cellOrder[nPathIdx+1];
-
-        // For every pair of consecutive cells, we search for a common face between them
-        for (label nFaceIdx : mesh.cells()[nCurrCellIdx])
-        {
-            label nNeiIdx = getNei(nCurrCellIdx, nFaceIdx);
-            // If there is a common face between them, we add a hit!
-            if (nNeiIdx == nNextCellIdx)
-            {
-                nHitCnt++;
-                break;
-            }
-        }
-    }
-    // The accuracy is the percentage of consecutive cells that were hits
-    return 100.0f * float(nHitCnt) / float(label(cellOrder.size()) - 1);
-}
-
-
-// Private methods
-void Foam::hpathRenumber::hpathFinder::initialize()
-{
-    // Data structure to keep track of cells already in the renumbering
-    bIsRenumbered = Foam::bitSet(mesh.nCells(), false);
-
-    // List used to seperate cells into layers
-    // Saves for every cell its layer index
-    nCellLayerIndex = Foam::labelList(mesh.nCells(), -1);
-
-    // List used to seperate cells within the same layer into seperate connected components
-    // Saves for every cell its connected component index
-    nCellConnnectedComponentIndex = Foam::labelList(mesh.nCells(), -1);
-
-    // Marks cells that have already been by BFS
-    bBFSFoundCell = Foam::bitSet(mesh.nCells(), false);
-    // Saves for every cell its face-neighbours within the same connected component
-    nConnnectedComponentGraph = Foam::labelListList(mesh.nCells());
-
-    // Saves for every cell within a layer its point-distance from the starting cell of that layer
-    nCellPointDistFromStart = Foam::labelList(mesh.nCells(), -1);
-    // Saves for every cell within a layer its face-distance from the starting cell of that layer
-    nCellFaceDistFromStart = Foam::labelList(mesh.nCells(), -1);
-
-    // For each cell its DFS depth within the connected component
-    nDFSCellDepth = Foam::labelList(mesh.nCells(), -1);
-    // For each cell its DFS parent within the connected component
-    nDFSParentCell = Foam::labelList(mesh.nCells(), -1);
-
-    // Now we can find the Mesh-Graph
-    nMeshGraph = Foam::labelListList(mesh.nCells());
-}
-
-
-Foam::label Foam::hpathRenumber::hpathFinder::getNei(label nCell, label nFaceIdx) const
-{
-    // If the face has no neighbor, return -1
-    if (nFaceIdx >= mesh.faceNeighbour().size())
-    {
-        return -1;
-    }
-
-    // Otherwise, the face connects 2 cells: the owner and the neighbor.
-    // One of these should be nCell
-    label nOwner = mesh.faceOwner()[nFaceIdx];
-    label nNei = mesh.faceNeighbour()[nFaceIdx];
-
-    // Find which of these two cells is the input cell, return the other one
-    return (nOwner == nCell) ? nNei : nOwner;
-}
-
-
-void Foam::hpathRenumber::hpathFinder::getMeshGraph()
-{
-    // The purpose of the bitSet is to avoid adding the same cell as a neighbour more than once
-    Foam::bitSet foundCells(mesh.nCells(), false);
-
-    for (label nCellIdx = 0; nCellIdx < mesh.nCells(); nCellIdx++)
-    {
-        // Dynamic list of point-neighbours of the current cell
-        DynamicList<label> nNeiList;
-
-        for (label nPntIdx : mesh.cellPoints()[nCellIdx])
-        {
-            for (label nNeiCell : mesh.pointCells()[nPntIdx])
-            {
-                if (nNeiCell == nCellIdx) continue;
-                if (foundCells[nNeiCell]) continue;
-
-                // We have found a cell with a common point to the current cell
-                // Therefore, we can now add it as a neighbour in the Mesh-Graph
-                nNeiList.append(nNeiCell);
-
-                foundCells[nNeiCell] = true;
-            }
-        }
-
-        // Convert from DynamicList to labelList
-        nMeshGraph[nCellIdx] = nNeiList;
-        
-        // Clear the bitSet before next iteration
-        for (label neiCell : nMeshGraph[nCellIdx])
-        {
-            foundCells[neiCell] = false;
-        }
-    }
-}
-
-
-void Foam::hpathRenumber::hpathFinder::getLayerSeparation
-(
-    Foam::DynamicList<Foam::DynamicList<label>>& nCellsByLayer
-)
-{
-    // We want to separate the mesh into layers:
-    //      1) The mesh is separated into connected components                                          - getConnectedComponents()
-    //      2) For each connected component we find the deepest cell and set it to be a starting cell   - getStartingCells()
-    //      3) Each connected component is separated into layers using a BFS from the starting          - getLayers()
-
-    std::cout << "Finding Layer Separation" << std::endl;
-
-    // Step 1: Finding connected components
-    Foam::labelList nAllCells(Foam::identity(mesh.nCells()));
-    Foam::DynamicList<Foam::DynamicList<label>> nCellsByConnnectedComponent;
-    getConnectedComponents(nAllCells, nCellsByConnnectedComponent);
-
-    // Step 2: Finding starting cells
-    // Finds for each connected component its deepest cell and returns them
-    // The starting cells are returned through nStartingCellList
-    // There will be exactly one starting cell per connected component
-    Foam::DynamicList<label> nStartingCellList;
-    getStartingCells(nStartingCellList);
-
-    // Step 3: Layer separation
-    // Layer separation is done in each connected component from the starting cell outwards
-    getLayers(nStartingCellList, nCellsByLayer);
-
-    label nLayerCnt = nCellsByLayer.size();
-    std::cout << "Mesh Separated into " << nLayerCnt << " layers" << std::endl;
-
-    // Reset connected component index list for solveLayer() to use
-    nCellConnnectedComponentIndex = Foam::labelList(mesh.nCells(), -1);
-}
-
-
-void Foam::hpathRenumber::hpathFinder::getStartingCells
-(
-    Foam::DynamicList<label>& nStartingCells)
-const
-{
-    // The starting cell in each connected component should be the 'deepest' cell in that connected component
-    // A cells 'depth' is its point-distance from any boundary cell
-
-    // The Algorithm:
-    //      1) Find all boundary points
-    //      2) Find all boundary cells
-    //      3) Separate the boundary cells by connected component
-    //      4) Find for each cell in the mesh its point-distance from the boundary
-    //      5) Return through nStartingCells a list containing the deepest cell from each connected component
-
-    // nCellDepthList will contain for each cell its depth within its ConnectedComponent
-    Foam::labelList nCellDepthList(mesh.nCells(), -1);
-
-    // Step 1: Finding all boundary points
-    //      - A boundary point is a point on a boundary face
-    Foam::bitSet bIsBoundaryPts(mesh.nPoints(), false);
-
-    // Iterate over boundary faces and mark their points as boundary
-    //  - Ignore boundary faces of type "empty"
-    const polyBoundaryMesh& bndMesh = mesh.boundaryMesh();
-    for (label nBndType = 0; nBndType < bndMesh.size(); ++nBndType)
-    {
-        // If boundary is of type "empty" - skip it
-        if (bndMesh[nBndType].type().compare("empty") == 0) continue;
-
-        labelRange range = bndMesh.patchRanges()[nBndType];
-        // For every face in range set all of its points as boundary
-        for (label nFaceIdx = range.min(); nFaceIdx <= range.max(); ++nFaceIdx)
-        {
-            for (label nPointIdx : mesh.faces()[nFaceIdx])
-            {
-                bIsBoundaryPts[nPointIdx] = true;
-            }
-        }
-    }
-    // Step 2: Finding all boundary cells
-    //      - A boundary cell is a cell containing at least one boundary point
-    Foam::bitSet bIsBoundaryCells(mesh.nCells(), false);
-
-    // For every boundary point: mark all cells that have it as boundary cells
-    for (label nPntIdx = 0; nPntIdx < mesh.nPoints(); nPntIdx++)
-    {
-        if (!bIsBoundaryPts[nPntIdx]) continue;
-        // Here we make use of the nPntCellList we found when computing the Mesh-Graph
-        for (label nCellIdx : mesh.pointCells()[nPntIdx])
-        {
-            bIsBoundaryCells[nCellIdx] = true;
-        }
-    }
-
-    // Step 3: Separate the boundary cells based on which connected component they are a part of
-    //      - Find for each connected component a list of all its boundary cells
-    label nConnnectedComponentCnt = *std::max_element(nCellConnnectedComponentIndex.begin(), nCellConnnectedComponentIndex.end()) + 1;
-    Foam::List<Foam::DynamicList<label>> nBoundaryCellsByConnnectedComponent(nConnnectedComponentCnt);
-    for (label nCellIdx = 0; nCellIdx < mesh.nCells(); nCellIdx++)
-    {
-        if (bIsBoundaryCells[nCellIdx])
-        {
-            nBoundaryCellsByConnnectedComponent[nCellConnnectedComponentIndex[nCellIdx]].append(nCellIdx);
-        }
-    }
-
-    // Steps 4+5:
-    //      - Find each cells distance from the boundary by running a BFS algorithm from the boundary in each connected component
-    //      - Along the way, find the deepest cell in each connected component and push them to the list
-
-    for (Foam::labelList& nBndCells : nBoundaryCellsByConnnectedComponent)
-    {
-        // We want to find the maximum depth cell within the connected component
-        label nMaxDepthCell = nBndCells[0];
-
-        // Run a BFS algorithm from the boundary:
-        //      - All boundary cells are pushed to the queue with depth 0
-        Foam::CircularBuffer<label> nBfsQueue;
-        for (label nCellIdx : nBndCells)
-        {
-            nCellDepthList[nCellIdx] = 0;
-            nBfsQueue.push_back(nCellIdx);
-        }
-
-        while (!nBfsQueue.empty())
-        {
-            label nCurrCell = nBfsQueue.first();
-            nBfsQueue.pop_front();
-
-            // Check if current cell is the deeper than the maximum depth cell found so far. If it is, we update nMaxDepthCell
-            if (nCellDepthList[nCurrCell] > nCellDepthList[nMaxDepthCell])
-            {
-                nMaxDepthCell = nCurrCell;
-            }
-
-            // The BFS algorithm needs to be based on point-neghbours, meaning by using the Mesh-Graph
-            // Note that while cells in different connected components are never face-neighbours, but they may be point-neighbours (neighbours in the Mesh-Graph)
-            // Therefore, when pushing all point-neighbouring cells we need to check that they are in the same connected component
-            for (label nNeiCell : nMeshGraph[nCurrCell])
-            {
-                if (nCellDepthList[nNeiCell] != -1) continue;
-                if (nCellConnnectedComponentIndex[nNeiCell] != nCellConnnectedComponentIndex[nCurrCell]) continue;
-
-                nCellDepthList[nNeiCell] = nCellDepthList[nCurrCell] + 1;
-                nBfsQueue.push_back(nNeiCell);
-            }
-        }
-
-        // Finally, we can push the maximum depth cell we found
-        nStartingCells.append(nMaxDepthCell);
-    }
-}
-
-
-void Foam::hpathRenumber::hpathFinder::getLayers
-(
-    const Foam::labelList& nStartingCells,
-    Foam::DynamicList<Foam::DynamicList<label>>& nCellsByLayer
-)
-{
-    // Separates all cells in the mesh to layers
-    // In each connected component:
-    //      1) The starting cell is set as 'Layer 0'
-    //      2) All cells that are point-neighbours of the starting cell are set as layer 1
-    //      3) All cells that are point-neighbours of a cell in layer 1 are set as layer 2
-    //      4) Repeat step (3) with increasing layer indices until all cells have been found
-    // In this way, cells in each component are grouped in layers by their MINIMUM point-distance to their starting cell
-
-    // Note: In reality, the same BFS is run for all components at once. Because components are connected, this is equivalent
-
-    // Now we run BFS from starting cells
-    Foam::CircularBuffer<label> nBfsQueue;
-
-    // Push all starting cells to queue
-    for (label nCellIdx : nStartingCells)
-    {
-        nCellLayerIndex[nCellIdx] = 0;
-        nBfsQueue.push_back(nCellIdx);
-    }
-
-    // BFS algorithm will find for each cell its minimum distance in the Mesh-Graph from the corresponding starting cell
-    while (!nBfsQueue.empty())
-    {
-        label nCurrCell = nBfsQueue.first();
-        nBfsQueue.pop_front();
-
-        label nLayer = nCellLayerIndex[nCurrCell];
-        if (label(nCellsByLayer.size()) <= nCellLayerIndex[nCurrCell])
-        {
-            nCellsByLayer.append(Foam::labelList());
-        }
-        nCellsByLayer[nLayer].append(nCurrCell);
-
-        for (label nNeiCell : nMeshGraph[nCurrCell])
-        {
-            // If neighbour is in a different connected component, ignore it
-            if (nCellConnnectedComponentIndex[nNeiCell] != nCellConnnectedComponentIndex[nCurrCell]) continue;
-            // If neighbour hasn't been visited yet, set it's layer and push it:
-            if (nCellLayerIndex[nNeiCell] != -1) continue;
-            nCellLayerIndex[nNeiCell] = nCellLayerIndex[nCurrCell] + 1;
-            nBfsQueue.push_back(nNeiCell);
-        }
-    }
-}
-
-
-void Foam::hpathRenumber::hpathFinder::solveLayer
-(
-    Foam::labelList nCellsInLayer,
-    Foam::labelList& cellOrder
-)
-{
-    // Note: this method assumes all cells in the nCellsInLayer have the same 'layer index' in nCellLayerIndex
-
-    // General rundown of the algorithm:
-    //      1) Cells are separated into connected components
-    //      2) Each connected component is solved separately and its path is added to the renumbering
-    //      3) If there are still cells that have not been found, return to step (1)
-
-    while (!nCellsInLayer.empty())
-    {
-        // Step 1: separate the cells into connected components
-        //      - Connected components need to be connected by FACES (not points)
-        Foam::DynamicList<Foam::DynamicList<label>> nCellsByConnnectedComponent;
-        getConnectedComponents(nCellsInLayer, nCellsByConnnectedComponent);
-
-        label nOrigFoundCellCount = nFoundCellCount;
-
-        // Step 2: For each connected component we call getHpathinConnnectedComponent()
-        for (const Foam::labelList& nCellsInConnectedComponent : nCellsByConnnectedComponent)
-        {
-            solveConnectedComponent(nCellsInConnectedComponent, cellOrder);
-        }
-        // Find how many cells were still not found
-        label nRemainingCellCount = nCellsInLayer.size() - (nFoundCellCount - nOrigFoundCellCount);
-        if (nRemainingCellCount == 0) break;
-
-        // Step 3: Find all the cells in the layer that were not found yet
-        Foam::labelList nRemainingCells(nRemainingCellCount);
-        label nIndex = 0;
-        for (label nCellIdx : nCellsInLayer)
-        {
-            if (!bIsRenumbered[nCellIdx])
-            {
-                nRemainingCells[nIndex++] = nCellIdx;
-            }
-        }
-
-        resetCells(nRemainingCells);
-        nCellsInLayer = nRemainingCells;
-    }
-}
-
-
-void Foam::hpathRenumber::hpathFinder::getConnectedComponents
-(
-    const Foam::labelList& nCellList,
-    Foam::DynamicList<Foam::DynamicList<label>>& nCellsByConnnectedComponent
-)
-{
-    // This function separates cells in a certain layer into connected components
-    // Each connected component will get a unique index
-    //      - nCellsByConnnectedComponent[i] will contain a list of all cells with connected component index 'i'
-
-    // We accomplish this using a generalized DFS algorithm:
-    //      While there are still cells that haven't been found:
-    //          1) Choose some cell in the layer that hasn't been found yet
-    //          2) Find all cells connected to that cell (using a dfs stack)
-    //          3) Set all of these cells as a new connected component
-
-    // Notice: This method does not use the Mesh-Graph!
-    //   - This is because the Mesh-Graph is POINT-connected, but we are searching for FACE-connected components
-
-    // The index of the current connected component
-    label nConnnectedComponentIndex = -1;
-
-    for (label nCellIdx : nCellList)
-    {
-        if (nCellConnnectedComponentIndex[nCellIdx] != -1) continue;
-
-        // This cell hasn't been found yet
-        // So we set it to be the beginning of a new connected component
-        // We increment nConnnectedComponentIndex, it is now the index of this new connected component
-        nConnnectedComponentIndex++;
-        nCellConnnectedComponentIndex[nCellIdx] = nConnnectedComponentIndex;
-        nCellsByConnnectedComponent.append(Foam::labelList());
-
-        Foam::CircularBuffer<label> nDfsStack;
-        nDfsStack.push_back(nCellIdx);
-
-        while(!nDfsStack.empty())
-        {
-            label nCurrCell = nDfsStack.last();
-            nDfsStack.pop_back();
-
-            nCellsByConnnectedComponent[nConnnectedComponentIndex].append(nCurrCell);
-
-            // We push all face-neighbouring cells that:
-            //      1) Are in the same layer
-            //      2) Haven't been found yet
-            for (label nFaceIdx : mesh.cells()[nCurrCell])
-            {
-                label nNeiCell = getNei(nCurrCell, nFaceIdx);
-                if (nNeiCell < 0) continue;
-                if (nCellLayerIndex[nNeiCell] != nCellLayerIndex[nCurrCell]) continue;
-                if (nCellConnnectedComponentIndex[nNeiCell] != -1) continue;
-
-                nCellConnnectedComponentIndex[nNeiCell] = nConnnectedComponentIndex;
-                nDfsStack.push_back(nNeiCell);
-            }
-        }
-    }
-}
-
-
-void Foam::hpathRenumber::hpathFinder::solveConnectedComponent
-(
-    const Foam::labelList& nCellsInConnectedComponent,
-    Foam::labelList& cellOrder
-)
-{
-    // For small cases, find path manually (1-2 cells)
-    if (nCellsInConnectedComponent.size() <= 2) {
-        for (label nCellIdx : nCellsInConnectedComponent) {
-            cellOrder[nFoundCellCount++] = nCellIdx;
-            bIsRenumbered[nCellIdx] = true;
-        }
-        return;
-    }
-    // Get a graph representing the connected component
-    getConnnectedComponentGraph(nCellsInConnectedComponent);
-
-    // Get the starting cell
-    label nStartCell = getStartingCellInConnnectedComponent(nCellsInConnectedComponent);
-
-    // Reorder each cells neighbours in the ConnnectedComponent-Graph in descending order by distance from start cell
-    reorderDistFromStart(nStartCell, nCellsInConnectedComponent);
-
-    // Get a path through the connected component, using the ConnnectedComponent-Graph
-    findPath(nStartCell, cellOrder);
-}
-
-
-void Foam::hpathRenumber::hpathFinder::getConnnectedComponentGraph
-(
-    const Foam::labelList& nCellsInConnectedComponent
-)
-{
-    // This method computes the ConnnectedComponent-Graph
-    // Cells are connected in the ConnnectedComponent-Graph if:
-    //      a) They are in the same layer
-    //      b) They are face-connected in the mesh (different from Mesh-Graph - there it was point-connected)
-    //          - note that this implies that they are also in the same connected component (hence the name)
-
-    for (label nCellIdx : nCellsInConnectedComponent)
-    {
-        Foam::DynamicList<label> nNeiList;
-
-        // For every face on the cell, find its neighbour through that face (if there is one)
-        // If that neighbour:
-        //      1) Hasn't been added to the reordering yet
-        //      2) Is in the same layer
-        //      3) Is in the same connected component
-        // Then we add it as a neighbour in the ConnectedComponent-Graph
-        for (label nFaceIdx : mesh.cells()[nCellIdx])
-        {
-            label nNeiCell = getNei(nCellIdx, nFaceIdx);
-            if (nNeiCell < 0) continue;
-            if (bIsRenumbered[nNeiCell]) continue;
-            if (nCellLayerIndex[nNeiCell] != nCellLayerIndex[nCellIdx]) continue;
-
-            nNeiList.append(nNeiCell);
-        }
-
-        nConnnectedComponentGraph[nCellIdx] = nNeiList;
-    }
-}
-
-
-Foam::label Foam::hpathRenumber::hpathFinder::getStartingCellInConnnectedComponent
-(
-    const Foam::labelList& nCellsInConnectedComponent
-)
-{
-
-    // Strategy for choosing a starting cell in the connected component:
-    //      1) Start from an arbitrary cell within the connected component
-    //      2) Find the furthest cell from it: choose it as the starting cell
-
-    // The reasoning behind this strategy is as follows:
-    //      - In many cases the cells given may be of the general shape of a long straight path
-    //      - In these cases, by starting from the farthest cell we guarantee it will be at one of the edges of the path
-
-    // It does not matter which cell we start from, so arbitrarly start from first cell in the connected component
-    label nCellIdx = nCellsInConnectedComponent[0];
-
-    // Find farthest cell from nCellIdx
-    // This is done using a standard BFS algorithm
-    Foam::CircularBuffer<label> nBfsQueue;
-
-    bBFSFoundCell[nCellIdx] = true;
-    nBfsQueue.push_back(nCellIdx);
-
-    label nCurrCell = -1;
-
-    while(!nBfsQueue.empty())
-    {
-        nCurrCell = nBfsQueue.first();
-        nBfsQueue.pop_front();
-
-        // Push all cells that have not yet been found by the BFS
-        for (label nNeiCell : nConnnectedComponentGraph[nCurrCell])
-        {
-            if (bBFSFoundCell[nNeiCell]) continue;
-            bBFSFoundCell[nNeiCell] = true;
-            nBfsQueue.push_back(nNeiCell);
-        }
-    }
-
-    // nCurrCell should now be the farthest cell from the starting cell in the connected component
-    return nCurrCell;
-}
-
-
-void Foam::hpathRenumber::hpathFinder::reorderDistFromStart
-(
-    label nStartCell,
-    const Foam::labelList& nCellsInConnectedComponent
-)
-{
-    // findPath() tends to find much better results when each cell's neighbours are ordered in a specific way based on their face/point-distance to the starting cell
-
-    // First, this method finds for each cell in the connected component its face-distance from the start cell
-    // Secondly, this method finds for each cell in the connected component its point-distance from the start cell
-    // Finally, this method reorders each cells neighbours in the ConnnectedComponent-Graph accordingly
-
-    //  - Both of these are done using a BFS algorithm: the first using face-neigbours and the second using point-neighbours
-
-    Foam::CircularBuffer<label> nBfsQueue;
-
-    // Find the face-distance of all cells in the connected component from the starting cell
-    nCellFaceDistFromStart[nStartCell] = 0;
-    // Push starting cell to queue
-    nBfsQueue.push_back(nStartCell);
-
-    while (!nBfsQueue.empty())
-    {
-        label nCurrCell = nBfsQueue.first();
-        nBfsQueue.pop_front();
-
-        // Push all face-neighbours that haven't already had their face distance found
-        // Face-neighbours are saved in the ConnnectedComponent-Graph
-        for (label nNeiCell : nConnnectedComponentGraph[nCurrCell])
-        {
-            if (nCellFaceDistFromStart[nNeiCell] != -1) continue;
-
-            nCellFaceDistFromStart[nNeiCell] = nCellFaceDistFromStart[nCurrCell] + 1;
-            nBfsQueue.push_back(nNeiCell);
-        }
-    }
-    // Find the *point*-distance of all cells in the connected component from the starting cell
-    nCellPointDistFromStart[nStartCell] = 0;
-    // Push starting cell to queue
-    nBfsQueue.push_back(nStartCell);
-
-    while (!nBfsQueue.empty())
-    {
-        label nCurrCell = nBfsQueue.first();
-        nBfsQueue.pop_front();
-
-        // Push all the point-neighbour that:
-        //      1) Haven't already been pushed by the BFS previously
-        //      2) Are in the same layer as the Starting Cell
-        //      3) Are in the same connected component as the Starting Cell
-        // Point neighbours are saved in the Mesh-Graph
-        for (label nNeiCell : nMeshGraph[nCurrCell])
-        {
-            if (bIsRenumbered[nNeiCell]) continue;
-            if (nCellPointDistFromStart[nNeiCell] != -1) continue;
-            if (nCellLayerIndex[nNeiCell] != nCellLayerIndex[nStartCell]) continue;
-            if (nCellConnnectedComponentIndex[nNeiCell] != nCellConnnectedComponentIndex[nStartCell]) continue;
-
-            nCellPointDistFromStart[nNeiCell] = nCellPointDistFromStart[nCurrCell] + 1;
-            nBfsQueue.push_back(nNeiCell);
-        }
-    }
-
-    for (label nCellIdx : nCellsInConnectedComponent)
-    {
-        // Sorting of neighbours is done in two levels:
-        //      1) Neigbours are sorted *descending*-order based on their *point*-distance from the starting cell
-        //      2) Neigbours with the same *point*-distance are sorted in *ascending*-order based on their *face*-distance from the starting cell
-        std::sort
-        (
-            nConnnectedComponentGraph[nCellIdx].begin(),
-            nConnnectedComponentGraph[nCellIdx].end(),
-            [this](label i, label j) {
-            if (nCellPointDistFromStart[i] != nCellPointDistFromStart[j])
-                return nCellPointDistFromStart[i] > nCellPointDistFromStart[j];
-            else
-                return nCellFaceDistFromStart[i] < nCellFaceDistFromStart[j];
-            }
-        );
-    }
-}
-
-
-void Foam::hpathRenumber::hpathFinder::findPath
-(
-    label nStartCell, Foam::labelList& cellOrder
-)
-{
-    // Tries to find the furthest cell from the starting cell within the ConnectedComponent-Graph
-    // This is done using a DFS algorithm:
-    //      - Run DFS from starting cell
-    //      - Return path to deepest cell found by DFS
-
-    // Run a standard DFS algorithm from starting cell
-    Foam::CircularBuffer<label> nDfsStack;
-    nDfsStack.push_back(nStartCell);
-    nDFSParentCell[nStartCell] = nStartCell;
-
-    // Used for finding and returning the best path
-    label nBestCell = nStartCell;
-
-    while (!nDfsStack.empty())
-    {
-        label nCurrCell = nDfsStack.last();
-        nDfsStack.pop_back();
-
-        // If the cell has already been popped previously, we can skip it
-        if (nDFSCellDepth[nCurrCell] != -1) continue;
-
-        // Otherwise, we update its depth value. This means the cell will never be PUSHED again
-        nDFSCellDepth[nCurrCell] = nDFSCellDepth[nDFSParentCell[nCurrCell]] + 1;
-
-        // If this is the new deepest cell, update best
-        if (nDFSCellDepth[nCurrCell] > nDFSCellDepth[nBestCell]) {
-            nBestCell = nCurrCell;
-        }
-
-        // Cells will be popped from the stack in reverse order from how we pushed them
-        // By iterating over the neighbors in reverse order, cells will be popped in the correct order
-        for (label nNeiIdx = nConnnectedComponentGraph[nCurrCell].size() - 1; nNeiIdx >= 0; nNeiIdx--)
-        {
-            label nNeiCell = nConnnectedComponentGraph[nCurrCell][nNeiIdx];
-            if (nDFSCellDepth[nNeiCell] == -1)
-            {
-                // Notice we only update the nCellDepth value of a cell when we POP it from the stack
-                // This means some cells may be pushed many times, but they will only be popped once
-                nDFSParentCell[nNeiCell] = nCurrCell;
-                nDfsStack.push_back(nNeiCell);
-            }
-        }
-    }
-
-    label nCellIdx = nBestCell;
-    label nPrevCell = -1;
-    while (nCellIdx != nPrevCell)
-    {
-        cellOrder[nFoundCellCount++] = nCellIdx;
-        bIsRenumbered[nCellIdx] = true;
-
-        // The first cell in the path is always its own parent
-        nPrevCell = nCellIdx;
-        nCellIdx = nDFSParentCell[nCellIdx];
-    }
-}
-
-
-void Foam::hpathRenumber::hpathFinder::resetCells
-(
-    Foam::labelList& nCellList
-)
-{
-    // Reset the data structures for the cells that were not found in the renumbering
-    for (label nCellIdx : nCellList)
-    {
-        nCellConnnectedComponentIndex[nCellIdx] = -1;
-        nCellPointDistFromStart[nCellIdx] = -1;
-        nCellFaceDistFromStart[nCellIdx] = -1;
-        bBFSFoundCell[nCellIdx] = false;
-        nDFSCellDepth[nCellIdx] = -1;
-        nDFSParentCell[nCellIdx] = -1;
-    }
-}
-
-
-// ************************************************************************* //

src/renumber/renumberMethods/hpathRenumber/hpathRenumber.H

@@ -1,142 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2023 Alon Zameret, Noam Manaker Morag
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::hpathRenumber
-
-Description
-    Renumbering based on Hamiltonian path.
-
-SourceFiles
-    hpathRenumber.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef Foam_hpathRenumber_H
-#define Foam_hpathRenumber_H
-
-#include "renumberMethod.H"
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                        Class hpathRenumber Declaration
-\*---------------------------------------------------------------------------*/
-
-class hpathRenumber
-:
-    public renumberMethod
-{
-    // Private Data
-
-        //- Apply layer separation? (default: true)
-        bool applyLayerSeparation_;
-
-
-    // Private Classes
-
-        // Helper class
-        class hpathFinder;
-
-public:
-
-    // Generated Methods
-
-        //- No copy construct
-        hpathRenumber(const hpathRenumber&) = delete;
-
-        //- No copy assignment
-        void operator=(const hpathRenumber&) = delete;
-
-
-    //- Runtime type information
-    TypeName("hpath");
-
-
-    // Constructors
-
-        //- Construct given the renumber dictionary
-        explicit hpathRenumber(const dictionary& dict);
-
-
-    //- Destructor
-    virtual ~hpathRenumber() = default;
-
-
-    // Member Functions
-
-        //- Return the order in which cells need to be visited
-        //- (ie. from ordered back to original cell label).
-        virtual labelList renumber(const pointField&) const
-        {
-            NotImplemented;
-            return labelList();
-        }
-
-        //- Return the order in which cells need to be visited
-        //- (ie. from ordered back to original cell label).
-        //  Uses mesh for regIOobject
-        virtual labelList renumber
-        (
-            const polyMesh& mesh,
-            const pointField& cc
-        ) const;
-
-        //- Return the order in which cells need to be visited
-        //- (ie. from ordered back to original cell label).
-        virtual labelList renumber
-        (
-            const CompactListList<label>& cellCells,
-            const pointField& cellCentres
-        ) const
-        {
-            NotImplemented;
-            return labelList();
-        }
-
-        //- Return the order in which cells need to be visited
-        //- (ie. from ordered back to original cell label).
-        virtual labelList renumber
-        (
-            const labelListList& cellCells,
-            const pointField& cellCentres
-        ) const
-        {
-            NotImplemented;
-            return labelList();
-        }
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/sampling/meshToMesh/meshToMeshTemplates.C

@@ -834,15 +834,15 @@ Foam::meshToMesh::mapTgtToSrc
         label srcPatchi = srcPatchID_[i];
         label tgtPatchi = tgtPatchID_[i];
 
-        if (!srcPatchFields.set(tgtPatchi))
+        if (!srcPatchFields.set(srcPatchi))
         {
             srcPatchFields.set
             (
                 srcPatchi,
                 fvPatchField<Type>::New
                 (
-                    tgtBfld[srcPatchi],
-                    srcMesh.boundary()[tgtPatchi],
+                    tgtBfld[tgtPatchi],
+                    srcMesh.boundary()[srcPatchi],
                     DimensionedField<Type, volMesh>::null(),
                     directFvPatchFieldMapper
                     (

src/surfMesh/readers/boundary/boundaryDataSurfaceReaderTemplates.C

@@ -88,7 +88,7 @@ Foam::boundaryDataSurfaceReader::readField
 {
     refPtr<Time> timePtr(Time::New(argList::envGlobalPath()));
 
-    return readField<Type>(baseDir, timeDir, fieldName, avg);
+    return readField<Type>(*timePtr, baseDir, timeDir, fieldName, avg);
 }
 
 

src/surfMesh/readers/ensight/ensightSurfaceReaderTemplates.C

@@ -82,7 +82,7 @@ Foam::tmp<Foam::Field<Type>> Foam::ensightSurfaceReader::readField
         }
 
         // Check that data type is as expected
-        // (assumes OpenFOAM generated the data set)
+        // (assuming OpenFOAM generated the data set)
         string primitiveType;
         is.read(primitiveType);
 
@@ -90,7 +90,8 @@ Foam::tmp<Foam::Field<Type>> Foam::ensightSurfaceReader::readField
 
         if
         (
-            primitiveType != ensightPTraits<Type>::typeName
+            debug
+         && primitiveType != ensightPTraits<Type>::typeName
          && primitiveType != pTraits<Type>::typeName
         )
         {

src/surfMesh/surfaceFormats/nas/NASsurfaceFormat.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2017-2022 OpenCFD Ltd.
+    Copyright (C) 2017-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -149,7 +149,6 @@ bool Foam::fileFormats::NASsurfaceFormat<Face>::read
     string line;
     while (is.good())
     {
-        string::size_type linei = 0;  // Parsing position within current line
         is.getLine(line);
 
         // ANSA extension
@@ -223,16 +222,30 @@ bool Foam::fileFormats::NASsurfaceFormat<Face>::read
             }
         }
 
+
+        // Parsing position within current line
+        std::string::size_type linei = 0;
+
+        // Is free format if line contains a comma
+        const bool freeFormat = line.contains(',');
+
         // First word (column 0-8)
         const word cmd(word::validate(nextNasField(line, linei, 8)));
 
         if (cmd == "CTRIA3")
         {
-            label elemId = readLabel(nextNasField(line, linei, 8)); // 8-16
-            label groupId = readLabel(nextNasField(line, linei, 8)); // 16-24
-            const auto a = readLabel(nextNasField(line, linei, 8)); // 24-32
-            const auto b = readLabel(nextNasField(line, linei, 8)); // 32-40
-            const auto c = readLabel(nextNasField(line, linei, 8)); // 40-48
+            // Fixed format:
+            //  8-16 : element id
+            // 16-24 : group id
+            // 24-32 : vertex
+            // 32-40 : vertex
+            // 40-48 : vertex
+
+            label elemId = readLabel(nextNasField(line, linei, 8, freeFormat));
+            label groupId = readLabel(nextNasField(line, linei, 8, freeFormat));
+            const auto a = readLabel(nextNasField(line, linei, 8, freeFormat));
+            const auto b = readLabel(nextNasField(line, linei, 8, freeFormat));
+            const auto c = readLabel(nextNasField(line, linei, 8, freeFormat));
 
             // Convert groupId into zoneId
             const auto iterZone = zoneLookup.cfind(groupId);
@@ -261,12 +274,20 @@ bool Foam::fileFormats::NASsurfaceFormat<Face>::read
         }
         else if (cmd == "CQUAD4")
         {
-            label elemId = readLabel(nextNasField(line, linei, 8)); // 8-16
-            label groupId = readLabel(nextNasField(line, linei, 8)); // 16-24
-            const auto a = readLabel(nextNasField(line, linei, 8)); // 24-32
-            const auto b = readLabel(nextNasField(line, linei, 8)); // 32-40
-            const auto c = readLabel(nextNasField(line, linei, 8)); // 40-48
-            const auto d = readLabel(nextNasField(line, linei, 8)); // 48-56
+            // Fixed format:
+            //  8-16 : element id
+            // 16-24 : group id
+            // 24-32 : vertex
+            // 32-40 : vertex
+            // 40-48 : vertex
+            // 48-56 : vertex
+
+            label elemId = readLabel(nextNasField(line, linei, 8, freeFormat));
+            label groupId = readLabel(nextNasField(line, linei, 8, freeFormat));
+            const auto a = readLabel(nextNasField(line, linei, 8, freeFormat));
+            const auto b = readLabel(nextNasField(line, linei, 8, freeFormat));
+            const auto c = readLabel(nextNasField(line, linei, 8, freeFormat));
+            const auto d = readLabel(nextNasField(line, linei, 8, freeFormat));
 
             // Convert groupId into zoneId
             const auto iterZone = zoneLookup.cfind(groupId);
@@ -310,11 +331,21 @@ bool Foam::fileFormats::NASsurfaceFormat<Face>::read
         }
         else if (cmd == "GRID")
         {
-            label index = readLabel(nextNasField(line, linei, 8)); // 8-16
-            (void) nextNasField(line, linei, 8); // 16-24
-            scalar x = readNasScalar(nextNasField(line, linei, 8)); // 24-32
-            scalar y = readNasScalar(nextNasField(line, linei, 8)); // 32-40
-            scalar z = readNasScalar(nextNasField(line, linei, 8)); // 40-48
+            // Fixed (short) format:
+            //  8-16 : point id
+            // 16-24 : coordinate system (not supported)
+            // 24-32 : point x coordinate
+            // 32-40 : point y coordinate
+            // 40-48 : point z coordinate
+            // 48-56 : displacement coordinate system (optional, unsupported)
+            // 56-64 : single point constraints (optional, unsupported)
+            // 64-70 : super-element id (optional, unsupported)
+
+            label index = readLabel(nextNasField(line, linei, 8, freeFormat));
+            (void) nextNasField(line, linei, 8, freeFormat);
+            scalar x = readNasScalar(nextNasField(line, linei, 8, freeFormat));
+            scalar y = readNasScalar(nextNasField(line, linei, 8, freeFormat));
+            scalar z = readNasScalar(nextNasField(line, linei, 8, freeFormat));
 
             pointId.append(index);
             dynPoints.append(point(x, y, z));
@@ -327,6 +358,8 @@ bool Foam::fileFormats::NASsurfaceFormat<Face>::read
             // GRID*      126   0 -5.55999875E+02 -5.68730474E+02
             // *         2.14897901E+02
 
+            // Cannot be long format and free format at the same time!
+
             label index = readLabel(nextNasField(line, linei, 16)); // 8-24
             (void) nextNasField(line, linei, 16); // 24-40
             scalar x = readNasScalar(nextNasField(line, linei, 16)); // 40-56

src/surfMesh/writers/nastran/nastranSurfaceWriter.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2022 OpenCFD Ltd.
+    Copyright (C) 2015-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -307,16 +307,10 @@ void Foam::surfaceWriters::nastranWriter::writeGeometry
 Foam::surfaceWriters::nastranWriter::nastranWriter()
 :
     surfaceWriter(),
-    writeFormat_(fieldFormat::SHORT),
-    fieldMap_(),
+    writeFormat_(fieldFormat::FREE),
     commonGeometry_(false),
-    separator_()
-{
-    // if (writeFormat_ == fieldFormat::FREE)
-    // {
-    //     separator_ = ",";
-    // }
-}
+    separator_(",")  // FREE format
+{}
 
 
 Foam::surfaceWriters::nastranWriter::nastranWriter
@@ -331,12 +325,10 @@ Foam::surfaceWriters::nastranWriter::nastranWriter
         (
             "format",
             options,
-            fieldFormat::LONG
+            fieldFormat::FREE
         )
     ),
-    fieldMap_(),
-    commonGeometry_(options.getOrDefault("commonGeometry", false)),
-    separator_()
+    commonGeometry_(options.getOrDefault("commonGeometry", false))
 {
     if (writeFormat_ == fieldFormat::FREE)
     {

src/surfMesh/writers/nastran/nastranSurfaceWriter.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2022 OpenCFD Ltd.
+    Copyright (C) 2015-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -33,13 +33,13 @@ Description
     The formatOptions for nastran:
     \table
         Property    | Description                           | Reqd | Default
-        fields      | Field pairs for PLOAD2/PLOAD4         | yes |
-        format      | Nastran format (short/long/free)      | no  | long
+        format      | Nastran format (short/long/free)      | no  | free
         scale       | Output geometry scaling               | no  | 1
         transform   | Output coordinate transform           | no  |
         fieldLevel  | Subtract field level before scaling   | no  | empty dict
         fieldScale  | Output field scaling                  | no  | empty dict
         commonGeometry | use separate geometry files        | no  | false
+        fields      | Field pairs for PLOAD2/PLOAD4         | yes |
     \endtable
 
     For example,
@@ -48,13 +48,6 @@ Description
     {
         nastran
         {
-            // OpenFOAM field name to NASTRAN load types
-            fields
-            (
-                (pMean  PLOAD2)
-                (p      PLOAD4)
-            );
-
             format  free;     // format type
 
             scale   1000;     // [m] -> [mm]
@@ -62,6 +55,13 @@ Description
             {
                "p.*"   0.01;  // [Pa] -> [mbar]
             }
+
+            // OpenFOAM field name to NASTRAN load types
+            fields
+            (
+                (pMean  PLOAD2)
+                (p      PLOAD4)
+            );
         }
     }
     \endverbatim
@@ -93,7 +93,6 @@ Description
 
 Note
     Output variable scaling does not apply to integer types such as Ids.
-    Field pairs default to PLOAD2 for scalars and PLOAD4 for vectors etc.
 
 SourceFiles
     nastranSurfaceWriter.C
@@ -221,10 +220,10 @@ public:
 
     // Constructors
 
-        //- Default construct. Default SHORT format
+        //- Default construct. Default FREE format
         nastranWriter();
 
-        //- Construct with some output options. Default LONG format
+        //- Construct with some output options. Default FREE format
         explicit nastranWriter(const dictionary& options);
 
         //- Construct from components

src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.C

@@ -128,7 +128,7 @@ Foam::scalar Foam::ReversibleReaction
     const scalarField& c
 ) const
 {
-    return kfwd/max(this->Kc(p, T), 1e-6);
+    return kfwd/max(this->Kc(p, T), VSMALL);
 }
 
 

src/transportModels/incompressible/viscosityModels/Casson/Casson.C

@@ -119,8 +119,8 @@ bool Foam::viscosityModels::Casson::read
 
     CassonCoeffs_.readEntry("m", m_);
     CassonCoeffs_.readEntry("tau0", tau0_);
-    CassonCoeffs_.readEntry("nuMin_", nuMin_);
-    CassonCoeffs_.readEntry("nuMax_", nuMax_);
+    CassonCoeffs_.readEntry("nuMin", nuMin_);
+    CassonCoeffs_.readEntry("nuMax", nuMax_);
 
     return true;
 }

tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/scalarTransport

@@ -7,7 +7,7 @@ tracer0
     log   off;
 
     resetOnStartUp  false;
-    // writeControl    writeTime;
+    writeControl    writeTime;
     // writeInterval   1;
     field           tracer0;
     D               0.001;

tutorials/heatTransfer/chtMultiRegionFoam/windshieldCondensation/system/controlDict

@@ -61,6 +61,8 @@ functions
         fvOptions
         {
         }
+
+        writeControl    writeTime;
     }
 
     fileUpdate

wmake/rules/General/cgal

@@ -1,21 +1,26 @@
 # ----------------------------------------------------------------------------
 # CGAL definitions - several possibilities
 #
-#  0. missing
-#  1. header-only
-#  2. library, no mpfr
-#  3. library, with mpfr (a likely default)
+#  - missing
+#  - header-only
+#  - header-only, no mpfr
+#  - library, no mpfr
+#  - library, with mpfr (default for older CGAL)
 #
 # Dispatch according to the defined 'CGAL_FLAVOUR'
 # - names may change [see wmake/scripts/have_cgal]
-# (no-cgal | cgal-header | cgal-no-mpfr | cgal-mpfr)
+# (no-cgal | cgal-header | cgal-header-no-mpfr | cgal-no-mpfr | cgal-mpfr)
 
 cgal_subrule := cgal-mpfr
-ifneq (,$(findstring header,$(CGAL_FLAVOUR)))
-    cgal_subrule := cgal-header-only
-endif
 ifneq (,$(findstring no-mpfr,$(CGAL_FLAVOUR)))
     cgal_subrule := cgal-no-mpfr
+    ifneq (,$(findstring header,$(CGAL_FLAVOUR)))
+        cgal_subrule := cgal-header-no-mpfr
+    endif
+else
+    ifneq (,$(findstring header,$(CGAL_FLAVOUR)))
+        cgal_subrule := cgal-header-only
+    endif
 endif
 
 # ----------------------------------------------------------------------------

wmake/rules/General/cgal-header-no-mpfr

@@ -0,0 +1,18 @@
+# -----------------------------------------------------------------------------
+# CGAL (header-only version) without mpfr
+
+CGAL_INC  = -DCGAL_HEADER_ONLY
+CGAL_LIBS =
+
+CGAL_INC  += \
+    $(foreach dir,$(BOOST_INC_DIR),-I$(dir)) \
+    $(foreach dir,$(CGAL_INC_DIR),-I$(dir))
+
+CGAL_LIBS += \
+    $(foreach dir,$(BOOST_LIB_DIR),-L$(dir))
+
+# ----
+# Extra failsafe - still needed? (2020-05-15)
+## CGAL_INC  += -I/usr/local/include -I/usr/include
+
+# -----------------------------------------------------------------------------

wmake/rules/darwin64Clang/cgal

@@ -1,14 +0,0 @@
-# ----------------------------------------------------------------------------
-# CGAL on Darwin
-# CGAL (library version) without mpfr
-
-CGAL_INC  = \
-    $(foreach dir,$(BOOST_INC_DIR),-I$(dir)) \
-    $(foreach dir,$(CGAL_INC_DIR),-I$(dir))
-
-CGAL_LIBS = \
-    $(foreach dir,$(BOOST_LIB_DIR),-L$(dir)) \
-    $(foreach dir,$(CGAL_LIB_DIR),-L$(dir)) \
-    -lCGAL
-
-# ----------------------------------------------------------------------------

wmake/rules/linux64Mingw/c++

@@ -35,6 +35,7 @@ LINKLIBSO   = $(CC) $(c++FLAGS) -shared \
     -Wl,--strip-all
 
 LINKEXE     = $(CC) $(c++FLAGS) \
+    -static-libgcc -static-libstdc++ \
     -Wl,--enable-auto-import \
     -Wl,--strip-all \
     -Wl,--force-exe-suffix