Changed code data structures to correspond with openfoam objects

ENH: new renumbering method based on Hamiltonian path
2023-02-27 12:07:41 +02:00 · 2023-02-27 10:04:47 +02:00 · 2023-01-31 16:25:57 +01:00 · 2023-01-17 13:51:52 +02:00
60 changed files with 1319 additions and 448 deletions
--- a/META-INFO/api-info
+++ b/META-INFO/api-info
@ -1,2 +1,2 @@
 api=2212
-patch=240625
+patch=230110
--- a/applications/solvers/compressible/rhoPimpleFoam/Make/options
+++ b/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 \
--- a/applications/solvers/compressible/rhoSimpleFoam/Make/options
+++ b/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 \
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options
+++ b/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 \
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/Make/options
+++ b/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 \
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options
+++ b/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 = \
--- a/applications/solvers/incompressible/pimpleFoam/Make/options
+++ b/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 \
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/CGALTriangulation3DKernel.H
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/CGALTriangulation3DKernel.H
@ -30,8 +30,8 @@ Description

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

-#ifndef Foam_CGALTriangulation3DKernel_H
-#define Foam_CGALTriangulation3DKernel_H
+#ifndef CGALTriangulation3DKernel_H
+#define CGALTriangulation3DKernel_H

 // Silence boost bind deprecation warnings (before CGAL-5.2.1)
 #include "CGAL/version.h"
@ -54,19 +54,9 @@ 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
--- a/applications/utilities/preProcessing/viewFactorsGen/Make/options
+++ b/applications/utilities/preProcessing/viewFactorsGen/Make/options
@ -4,6 +4,7 @@ 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 \
@ -11,6 +12,7 @@ EXE_INC = \


 EXE_LIBS = \
+    /* ${CGAL_LIBS} */ \
    -lfiniteVolume \
    -lsurfMesh \
    -lmeshTools \
--- a/etc/caseDicts/postProcessing/solvers/scalarTransport/scalarTransport
+++ b/etc/caseDicts/postProcessing/solvers/scalarTransport/scalarTransport
@ -18,8 +18,6 @@ Description
 type            scalarTransport;
 libs            ("libsolverFunctionObjects.so");

-writeControl    writeTime;
-
 field           s;
 schemesField    s;
 D               1e-09;
--- a/src/OpenFOAM/db/dictionary/functionEntries/evalEntry/evalEntry.C
+++ b/src/OpenFOAM/db/dictionary/functionEntries/evalEntry/evalEntry.C
@ -174,7 +174,7 @@ Foam::tokenList Foam::functionEntries::evalEntry::evaluate
        result.writeField(toks);
    }

-    return tokenList(std::move(toks.tokens()));
+    return std::move(toks);
 }


--- a/src/OpenFOAM/expressions/exprDriver/exprDriverFields.C
+++ b/src/OpenFOAM/expressions/exprDriver/exprDriverFields.C
@ -5,7 +5,7 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2023 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -41,11 +41,9 @@ void Foam::expressions::exprDriver::fill_random
 {
    if (seed <= 0)
    {
-        const TimeState* ts = this->timeState();
-
-        if (ts)
+        if (timeStatePtr_)
        {
-            seed = (ts->timeIndex() - seed);
+            seed = (timeStatePtr_->timeIndex() - seed);
        }
        else
        {
--- a/src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.C
+++ b/src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.C
@ -152,11 +152,10 @@ slicedBoundaryField
                new SlicedPatchField<Type>
                (
                    mesh.boundary()[patchi],
-                    DimensionedField<Type, GeoMesh>::null()
+                    DimensionedField<Type, GeoMesh>::null(),
+                    bField[patchi]
                )
            );
-
-            bf[patchi].UList<Type>::shallowCopy(bField[patchi]);
        }
    }

--- a/src/OpenFOAM/matrices/Matrix/Matrix.C
+++ b/src/OpenFOAM/matrices/Matrix/Matrix.C
@ -201,8 +201,7 @@ inline Foam::Matrix<Form, Type>::Matrix
 )
 :
    mRows_(Mb.m()),
-    nCols_(Mb.n()),
-    v_(nullptr)
+    nCols_(Mb.n())
 {
    doAlloc();

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

--- a/src/OpenFOAM/parallel/globalIndex/globalIndex.H
+++ b/src/OpenFOAM/parallel/globalIndex/globalIndex.H
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2023 OpenCFD Ltd.
+    Copyright (C) 2018-2022 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 ProcIDsContainer, class Type, class Addr>
+        template<class Type, class Addr>
        static void gather
        (
            const labelUList& offsets,  //!< offsets (master only)
            const label comm,           //!< communicator
-            const ProcIDsContainer& procIDs,
+            const UList<int>& procIDs,
            const IndirectListBase<Type, Addr>& fld,
            List<Type>& allFld,         //! output field (master only)
            const int tag = UPstream::msgType(),
--- a/src/OpenFOAM/parallel/globalIndex/globalIndexTemplates.C
+++ b/src/OpenFOAM/parallel/globalIndex/globalIndexTemplates.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2013-2017 OpenFOAM Foundation
-    Copyright (C) 2019-2023 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -248,12 +248,12 @@ void Foam::globalIndex::gather
 }


-template<class ProcIDsContainer, class Type, class Addr>
+template<class Type, class Addr>
 void Foam::globalIndex::gather
 (
    const labelUList& off,  // needed on master only
    const label comm,
-    const ProcIDsContainer& procIDs,
+    const UList<int>& 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::allProcs(comm),   // All communicator ranks
+            UPstream::procID(comm),
            sendData,
            allData,
            tag,
@ -404,7 +404,7 @@ void Foam::globalIndex::gather
        (
            offsets_,  // needed on master only
            comm,
-            UPstream::allProcs(comm),   // All communicator ranks
+            UPstream::procID(comm),
            sendData,
            allData,
            tag,
@ -622,7 +622,7 @@ void Foam::globalIndex::mpiGather
            (
                offsets_,  // needed on master only
                comm,
-                UPstream::allProcs(comm),   // All communicator ranks
+                UPstream::procID(comm),
                sendData,
                allData,
                tag,
@ -967,7 +967,7 @@ void Foam::globalIndex::scatter
        (
            offsets_,  // needed on master only
            comm,
-            UPstream::allProcs(comm),  // All communicator ranks
+            UPstream::procID(comm),
            allData,
            localData,
            tag,
--- a/src/OpenFOAM/primitives/Tensor/TensorI.H
+++ b/src/OpenFOAM/primitives/Tensor/TensorI.H
@ -409,10 +409,6 @@ 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
 {
@ -436,10 +432,6 @@ 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
 {
@ -875,10 +867,6 @@ 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)
 {
@ -893,10 +881,6 @@ 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)
 {
@ -911,10 +895,6 @@ 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)
 {
@ -937,10 +917,6 @@ 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)
 {
@ -981,10 +957,6 @@ 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)
 {
--- a/src/OpenFOAM/primitives/strings/word/word.C
+++ b/src/OpenFOAM/primitives/strings/word/word.C
@ -29,7 +29,6 @@ License
 #include "word.H"
 #include "debug.H"
 #include <cctype>
-#include <cstdint>
 #include <sstream>

 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
--- a/src/Pstream/mpi/UPstream.C
+++ b/src/Pstream/mpi/UPstream.C
@ -641,16 +641,7 @@ void Foam::UPstream::freePstreamCommunicator(const label communicator)
    }

    // Not touching the first two communicators (SELF, WORLD)
-    // 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 (communicator > 1)
    {
        if (MPI_COMM_NULL != PstreamGlobals::MPICommunicators_[communicator])
        {
--- a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/turbulentFluidThermoModels.C
+++ b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/turbulentFluidThermoModels.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2013-2016 OpenFOAM Foundation
-    Copyright (C) 2022-2023 OpenCFD Ltd.
+    Copyright (C) 2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -131,9 +131,6 @@ makeLESModel(dynamicKEqn);
 #include "dynamicLagrangian.H"
 makeLESModel(dynamicLagrangian);

-#include "sigma.H"
-makeLESModel(sigma);
-
 #include "SpalartAllmarasDES.H"
 makeLESModel(SpalartAllmarasDES);

--- a/src/TurbulenceModels/turbulenceModels/laminar/generalizedNewtonian/generalizedNewtonianViscosityModels/Casson/Casson.C
+++ b/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;
 }
--- a/src/fileFormats/nastran/NASCore.C
+++ b/src/fileFormats/nastran/NASCore.C
@ -5,7 +5,7 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2024 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -122,8 +122,7 @@ std::string Foam::fileFormats::NASCore::nextNasField
 (
    const std::string& str,
    std::string::size_type& pos,
-    const std::string::size_type width,
-    const bool free_format
+    std::string::size_type len
 )
 {
    const auto beg = pos;
@ -131,23 +130,15 @@ std::string Foam::fileFormats::NASCore::nextNasField

    if (end == std::string::npos)
    {
-        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);
+        pos = beg + len;    // Continue after field width
    }
    else
    {
-        // Free format - continue after comma
-        pos = end + 1;
-        return str.substr(beg, (end - beg));
+        len = (end - beg);  // Efffective width
+        pos = end + 1;      // Continue after comma
    }
+
+    return str.substr(beg, len);
 }


@ -244,8 +235,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 y coordinate
-    // 6 X3   : point z coordinate
+    // 5 X2   : point x coordinate
+    // 6 X3   : point x coordinate
    // 7 CD   : coordinate system for displacements (blank)
    // 8 PS   : single point constraints            (blank)
    // 9 SEID : super-element ID
--- a/src/fileFormats/nastran/NASCore.H
+++ b/src/fileFormats/nastran/NASCore.H
@ -5,7 +5,7 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2024 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -48,6 +48,7 @@ SourceFiles

 namespace Foam
 {
+
 namespace fileFormats
 {

@ -73,18 +74,18 @@ public:
        //- Output load format
        enum loadFormat
        {
-            PLOAD2,     //!< Face load (eg, pressure)
-            PLOAD4      //!< Vertex load
+            PLOAD2,
+            PLOAD4
        };

-        //- Selection names for the NASTRAN load formats
+        //- Selection names for the NASTRAN file field formats
        static const Enum<loadFormat> loadFormatNames;


    // Constructors

        //- Default construct
-        NASCore() noexcept = default;
+        NASCore() = default;


    // Public Static Member Functions
@ -92,20 +93,18 @@ public:
        //- Extract numbers from things like "-2.358-8" (same as "-2.358e-8")
        static scalar readNasScalar(const std::string& str);

-        //- 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
+        //- 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.
        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,
-            //! 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
+            std::string::size_type len
        );


--- a/src/fileFormats/stl/STLCore.C
+++ b/src/fileFormats/stl/STLCore.C
@ -135,9 +135,8 @@ int Foam::fileFormats::STLCore::detectBinaryHeader

        bad =
            (
-                dataFileSize < STLHeaderSize
-             || nTris < (dataFileSize - STLHeaderSize)/50
-             || nTris > (dataFileSize - STLHeaderSize)/25
+                nTris < int(dataFileSize - STLHeaderSize)/50
+             || nTris > int(dataFileSize - STLHeaderSize)/25
            );
    }

@ -209,9 +208,8 @@ Foam::fileFormats::STLCore::readBinaryHeader

        bad =
            (
-                dataFileSize < STLHeaderSize
-             || nTris < (dataFileSize - STLHeaderSize)/50
-             || nTris > (dataFileSize - STLHeaderSize)/25
+                nTris < int(dataFileSize - STLHeaderSize)/50
+             || nTris > int(dataFileSize - STLHeaderSize)/25
            );
    }

--- a/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.C
+++ b/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2018-2023 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -180,26 +180,12 @@ Foam::faBoundaryMesh::faBoundaryMesh

 void Foam::faBoundaryMesh::calcGeometry()
 {
-    // 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.
-
+    // 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.
-
-    (void)mesh_.edgeAreaNormals();
    (void)mesh_.pointAreaNormals();

-    (void)mesh_.areaCentres();
-    (void)mesh_.faceAreaNormals();
-
-
    PstreamBuffers pBufs(Pstream::defaultCommsType);

    if
@ -787,15 +773,12 @@ bool Foam::faBoundaryMesh::checkDefinition(const bool report) const

 void Foam::faBoundaryMesh::movePoints(const pointField& p)
 {
-    // See comments in calcGeometry()
-
-    (void)mesh_.edgeAreaNormals();
+    // 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.
    (void)mesh_.pointAreaNormals();

-    (void)mesh_.areaCentres();
-    (void)mesh_.faceAreaNormals();
-
-
    PstreamBuffers pBufs(Pstream::defaultCommsType);

    if
--- a/src/finiteArea/faMesh/faMesh.C
+++ b/src/finiteArea/faMesh/faMesh.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2020-2023 OpenCFD Ltd.
+    Copyright (C) 2020-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -277,27 +277,6 @@ 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)
@ -317,7 +296,18 @@ bool Foam::faMesh::init(const bool doInit)
    // Calculate the geometry for the patches (transformation tensors etc.)
    boundary_.calcGeometry();

-    syncGeom();
+    // 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>();
+    }

    return false;
 }
@ -999,6 +989,7 @@ bool Foam::faMesh::movePoints()

    clearGeomNotAreas();

+    // To satisfy the motion interface for MeshObject, const cast is needed
    if (patchPtr_)
    {
        patchPtr_->movePoints(newPoints);
@ -1012,8 +1003,6 @@ bool Foam::faMesh::movePoints()

    // Note: Fluxes were dummy?

-    syncGeom();
-
    return true;
 }

--- a/src/finiteArea/faMesh/faMesh.H
+++ b/src/finiteArea/faMesh/faMesh.H
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2021-2023 OpenCFD Ltd.
+    Copyright (C) 2021-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -617,10 +617,6 @@ 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

--- a/src/finiteArea/faMesh/faMeshDemandDrivenData.C
+++ b/src/finiteArea/faMesh/faMeshDemandDrivenData.C
@ -898,12 +898,6 @@ void Foam::faMesh::calcFaceCentres() const
            }
        }
    }
-
-    // Parallel consistency, exchange on processor patches
-    if (UPstream::parRun())
-    {
-        centres.boundaryFieldRef().evaluateCoupled<processorFaPatch>();
-    }
 }


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


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


--- a/src/functionObjects/field/expressions/fvExpressionField.C
+++ b/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("store", true);
+    store_ = dict.getOrDefault("autowrite", true);

    // "dimensions" is optional
    dimensions_.clear();
--- a/src/functionObjects/utilities/ensightWrite/ensightWriteUpdate.C
+++ b/src/functionObjects/utilities/ensightWrite/ensightWriteUpdate.C
@ -103,7 +103,7 @@ bool Foam::functionObjects::ensightWrite::readSelection(const dictionary& dict)
    selectFields_.uniq();

    blockFields_.clear();
-    dict.readIfPresent("excludeFields", blockFields_);
+    dict.readIfPresent("blockFields", blockFields_);
    blockFields_.uniq();

    // Actions to define selection
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectedParticleDistributionInjection/InjectedParticleDistributionInjection.C
+++ b/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/6.0;
+            const scalar volume = sumPow3*mathematical::pi/16.0;
            sumVolume += volume;
            volumeFlowRate_[injectori] = volume/dTime;

--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectedParticleInjection/InjectedParticleInjection.C
+++ b/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/6.0;
+        scalar vol = pow3(diameter_[i])*mathematical::pi/16.0;
        volume[i] = vol;
        sumVolume += vol;
    }
--- a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C
+++ b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2019-2023 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -62,53 +62,6 @@ 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
 (
--- a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.H
+++ b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.H
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2021-2023 OpenCFD Ltd.
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -143,21 +143,6 @@ 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
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C
@ -1368,7 +1368,6 @@ Foam::label Foam::snappyRefineDriver::refinementInterfaceRefine
                                (
                                    face2i != facei
                                 && surfaceIndex[face2i] != -1
-                                 && cutter.faceLevel(face2i) > cLevel
                                )
                                {
                                    // Get outwards pointing normal
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/faceAreaWeightAMI/faceAreaWeightAMI.C
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/faceAreaWeightAMI/faceAreaWeightAMI.C
@ -269,8 +269,7 @@ bool Foam::faceAreaWeightAMI::setNextFaces
        return false;
    }

-    const auto& srcPatch = this->srcPatch();
-    const labelList& srcNbrFaces = srcPatch.faceFaces()[srcFacei];
+    const labelList& srcNbrFaces = this->srcPatch().faceFaces()[srcFacei];

    // Initialise tgtFacei
    tgtFacei = -1;
@ -361,7 +360,6 @@ bool Foam::faceAreaWeightAMI::setNextFaces
    {
        FatalErrorInFunction
            << "Unable to set target face for source face " << srcFacei
-            << " with centre: " << srcPatch.faceCentres()[srcFacei]
            << abort(FatalError);
    }

--- a/src/meshTools/edgeMesh/edgeFormats/nastran/NASedgeFormat.C
+++ b/src/meshTools/edgeMesh/edgeFormats/nastran/NASedgeFormat.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2017-2024 OpenCFD Ltd.
+    Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -62,17 +62,13 @@ 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())
+        if (line.empty() || line[0] == '$')
        {
-            continue;  // Ignore empty
-        }
-        else if (line[0] == '$')
-        {
-            // Ignore comment
-            continue;
+            continue; // Skip empty or comment
        }

        // Check if character 72 is continuation
@ -98,66 +94,38 @@ 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")
        {
-            // Fixed format:
-            //  8-16 : element id
-            // 16-24 : group id
-            // 24-32 : vertex
-            // 32-40 : vertex
+            // discard elementId (8-16)
+            (void) nextNasField(line, linei, 8); // 8-16
+            // discard groupId (16-24)
+            (void) nextNasField(line, linei, 8); // 16-24

-            // 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));
+            label a = readLabel(nextNasField(line, linei, 8)); // 24-32
+            label b = readLabel(nextNasField(line, linei, 8)); // 32-40

            dynEdges.append(edge(a,b));
        }
        else if (cmd == "PLOTEL")
        {
-            // Fixed format:
-            //  8-16 : element id
-            // 16-24 : vertex
-            // 24-32 : vertex
-            // 32-40 : vertex
-
            // discard elementId (8-16)
-            (void) nextNasField(line, linei, 8, freeFormat);
+            (void) nextNasField(line, linei, 8); // 8-16

-            label a = readLabel(nextNasField(line, linei, 8, freeFormat));
-            label b = readLabel(nextNasField(line, linei, 8, freeFormat));
+            label a = readLabel(nextNasField(line, linei, 8)); // 16-24
+            label b = readLabel(nextNasField(line, linei, 8)); // 24-32

            dynEdges.append(edge(a,b));
        }
        else if (cmd == "GRID")
        {
-            // 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));
+            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

            pointId.append(index);
            dynPoints.append(point(x, y, z));
@ -170,8 +138,6 @@ 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
--- a/src/meshTools/output/vtk/topoSet/foamVtkWritePointSet.C
+++ b/src/meshTools/output/vtk/topoSet/foamVtkWritePointSet.C
@ -155,7 +155,7 @@ bool Foam::vtk::writePointSet

    if (parallel)
    {
-        vtk::writeListParallel(format.ref(), mesh.points(), pointLabels);
+        vtk::writeListParallel(format(), mesh.points(), pointLabels);
    }
    else
    {
--- a/src/optimisation/adjointOptimisation/adjoint/ATCModel/ATCUaGradU/ATCUaGradU.C
+++ b/src/optimisation/adjointOptimisation/adjoint/ATCModel/ATCUaGradU/ATCUaGradU.C
@ -5,8 +5,8 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2023 PCOpt/NTUA
-    Copyright (C) 2013-2023 FOSS GP
+    Copyright (C) 2007-2021 PCOpt/NTUA
+    Copyright (C) 2013-2021 FOSS GP
    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@ -86,9 +86,10 @@ void ATCUaGradU::addATC(fvVectorMatrix& UaEqn)
    if (extraConvection_ > 0)
    {
        // Implicit part added to increase diagonal dominance
-        UaEqn += ATClimiter_*extraConvection_*fvm::div(-phi, Ua);
+        // Note: Maybe this needs to be multiplied with the ATClimiter ??
+        UaEqn += 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);
    }

--- a/src/optimisation/adjointOptimisation/adjoint/ATCModel/ATCstandard/ATCstandard.C
+++ b/src/optimisation/adjointOptimisation/adjoint/ATCModel/ATCstandard/ATCstandard.C
@ -5,8 +5,8 @@
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2023 PCOpt/NTUA
-    Copyright (C) 2013-2023 FOSS GP
+    Copyright (C) 2007-2021 PCOpt/NTUA
+    Copyright (C) 2013-2021 FOSS GP
    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@ -90,9 +90,10 @@ void ATCstandard::addATC(fvVectorMatrix& UaEqn)
    if (extraConvection_ > 0)
    {
        // Implicit part added to increase diagonal dominance
-        UaEqn += ATClimiter_*extraConvection_*fvm::div(-phi, Ua);
+        // Note: Maybe this needs to be multiplied with the ATClimiter ??
+        UaEqn += 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);
    }

--- a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
+++ b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
@ -1949,7 +1949,7 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
 //        //globalBorderTris.gather
 //        //(
 //        //    UPstream::worldComm,
-//        //    UPstream::allProcs(UPstream::worldComm),
+//        //    UPstream::procID(Pstream::worldComm),
 //        //    globalBorderCentres
 //        //);
 //        pointField globalBorderCentres(allCentres);
@ -1996,7 +1996,7 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
 //        //globalBorderTris.scatter
 //        //(
 //        //    UPstream::worldComm,
-//        //    UPstream::allProcs(UPstream::worldComm),
+//        //    UPstream::procID(Pstream::worldComm),
 //        //    isMasterPoint
 //        //);
 //        //boolList isMasterBorder(s.size(), false);
@ -2094,7 +2094,7 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
            globalTris().gather
            (
                UPstream::worldComm,
-                UPstream::allProcs(UPstream::worldComm),
+                UPstream::procID(Pstream::worldComm),
                allCentres
            );
        }
@ -2144,7 +2144,7 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
            globalTris().scatter
            (
                UPstream::worldComm,
-                UPstream::allProcs(UPstream::worldComm),
+                UPstream::procID(Pstream::worldComm),
                allDistribution,
                distribution
            );
--- a/src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/alphatPhaseChangeWallFunction/alphatPhaseChangeWallFunctionFvPatchScalarField.C
+++ b/src/phaseSystemModels/reactingEuler/multiphaseSystem/derivedFvPatchFields/alphatPhaseChangeWallFunction/alphatPhaseChangeWallFunctionFvPatchScalarField.C
@ -75,7 +75,7 @@ alphatPhaseChangeWallFunctionFvPatchScalarField

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

--- a/src/renumber/renumberMethods/Make/files
+++ b/src/renumber/renumberMethods/Make/files
@ -5,5 +5,6 @@ randomRenumber/randomRenumber.C
 springRenumber/springRenumber.C
 structuredRenumber/structuredRenumber.C
 structuredRenumber/OppositeFaceCellWaveBase.C
+hpathRenumber/hpathRenumber.C

 LIB = $(FOAM_LIBBIN)/librenumberMethods
--- a/src/renumber/renumberMethods/hpathRenumber/hpathRenumber.C
+++ b/src/renumber/renumberMethods/hpathRenumber/hpathRenumber.C
@ -0,0 +1,976 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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;
+    }
+}
+
+
+// ************************************************************************* //
--- a/src/renumber/renumberMethods/hpathRenumber/hpathRenumber.H
+++ b/src/renumber/renumberMethods/hpathRenumber/hpathRenumber.H
@ -0,0 +1,142 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+// ************************************************************************* //
--- a/src/sampling/meshToMesh/meshToMeshTemplates.C
+++ b/src/sampling/meshToMesh/meshToMeshTemplates.C
@ -834,15 +834,15 @@ Foam::meshToMesh::mapTgtToSrc
        label srcPatchi = srcPatchID_[i];
        label tgtPatchi = tgtPatchID_[i];

-        if (!srcPatchFields.set(srcPatchi))
+        if (!srcPatchFields.set(tgtPatchi))
        {
            srcPatchFields.set
            (
                srcPatchi,
                fvPatchField<Type>::New
                (
-                    tgtBfld[tgtPatchi],
-                    srcMesh.boundary()[srcPatchi],
+                    tgtBfld[srcPatchi],
+                    srcMesh.boundary()[tgtPatchi],
                    DimensionedField<Type, volMesh>::null(),
                    directFvPatchFieldMapper
                    (
--- a/src/surfMesh/readers/boundary/boundaryDataSurfaceReaderTemplates.C
+++ b/src/surfMesh/readers/boundary/boundaryDataSurfaceReaderTemplates.C
@ -88,7 +88,7 @@ Foam::boundaryDataSurfaceReader::readField
 {
    refPtr<Time> timePtr(Time::New(argList::envGlobalPath()));

-    return readField<Type>(*timePtr, baseDir, timeDir, fieldName, avg);
+    return readField<Type>(baseDir, timeDir, fieldName, avg);
 }


--- a/src/surfMesh/readers/ensight/ensightSurfaceReaderTemplates.C
+++ b/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
-        // (assuming OpenFOAM generated the data set)
+        // (assumes OpenFOAM generated the data set)
        string primitiveType;
        is.read(primitiveType);

@ -90,8 +90,7 @@ Foam::tmp<Foam::Field<Type>> Foam::ensightSurfaceReader::readField

        if
        (
-            debug
-         && primitiveType != ensightPTraits<Type>::typeName
+            primitiveType != ensightPTraits<Type>::typeName
         && primitiveType != pTraits<Type>::typeName
        )
        {
--- a/src/surfMesh/surfaceFormats/nas/NASsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/nas/NASsurfaceFormat.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2017-2024 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -149,6 +149,7 @@ 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
@ -222,30 +223,16 @@ 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")
        {
-            // 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));
+            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

            // Convert groupId into zoneId
            const auto iterZone = zoneLookup.cfind(groupId);
@ -274,20 +261,12 @@ bool Foam::fileFormats::NASsurfaceFormat<Face>::read
        }
        else if (cmd == "CQUAD4")
        {
-            // 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));
+            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

            // Convert groupId into zoneId
            const auto iterZone = zoneLookup.cfind(groupId);
@ -331,21 +310,11 @@ bool Foam::fileFormats::NASsurfaceFormat<Face>::read
        }
        else if (cmd == "GRID")
        {
-            // 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));
+            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

            pointId.append(index);
            dynPoints.append(point(x, y, z));
@ -358,8 +327,6 @@ 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
--- a/src/surfMesh/writers/nastran/nastranSurfaceWriter.C
+++ b/src/surfMesh/writers/nastran/nastranSurfaceWriter.C
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2024 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -307,10 +307,16 @@ void Foam::surfaceWriters::nastranWriter::writeGeometry
 Foam::surfaceWriters::nastranWriter::nastranWriter()
 :
    surfaceWriter(),
-    writeFormat_(fieldFormat::FREE),
+    writeFormat_(fieldFormat::SHORT),
+    fieldMap_(),
    commonGeometry_(false),
-    separator_(",")  // FREE format
-{}
+    separator_()
+{
+    // if (writeFormat_ == fieldFormat::FREE)
+    // {
+    //     separator_ = ",";
+    // }
+}


 Foam::surfaceWriters::nastranWriter::nastranWriter
@ -325,10 +331,12 @@ Foam::surfaceWriters::nastranWriter::nastranWriter
        (
            "format",
            options,
-            fieldFormat::FREE
+            fieldFormat::LONG
        )
    ),
-    commonGeometry_(options.getOrDefault("commonGeometry", false))
+    fieldMap_(),
+    commonGeometry_(options.getOrDefault("commonGeometry", false)),
+    separator_()
 {
    if (writeFormat_ == fieldFormat::FREE)
    {
--- a/src/surfMesh/writers/nastran/nastranSurfaceWriter.H
+++ b/src/surfMesh/writers/nastran/nastranSurfaceWriter.H
@ -6,7 +6,7 @@
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
    Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2024 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
    This file is part of OpenFOAM.
@ -33,13 +33,13 @@ Description
    The formatOptions for nastran:
    \table
        Property    | Description                           | Reqd | Default
-        format      | Nastran format (short/long/free)      | no  | free
+        fields      | Field pairs for PLOAD2/PLOAD4         | yes |
+        format      | Nastran format (short/long/free)      | no  | long
        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,6 +48,13 @@ Description
    {
        nastran
        {
+            // OpenFOAM field name to NASTRAN load types
+            fields
+            (
+                (pMean  PLOAD2)
+                (p      PLOAD4)
+            );
+
            format  free;     // format type

            scale   1000;     // [m] -> [mm]
@ -55,13 +62,6 @@ Description
            {
               "p.*"   0.01;  // [Pa] -> [mbar]
            }
-
-            // OpenFOAM field name to NASTRAN load types
-            fields
-            (
-                (pMean  PLOAD2)
-                (p      PLOAD4)
-            );
        }
    }
    \endverbatim
@ -93,6 +93,7 @@ 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
@ -220,10 +221,10 @@ public:

    // Constructors

-        //- Default construct. Default FREE format
+        //- Default construct. Default SHORT format
        nastranWriter();

-        //- Construct with some output options. Default FREE format
+        //- Construct with some output options. Default LONG format
        explicit nastranWriter(const dictionary& options);

        //- Construct from components
--- a/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.C
+++ b/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), VSMALL);
+    return kfwd/max(this->Kc(p, T), 1e-6);
 }


--- a/src/transportModels/incompressible/viscosityModels/Casson/Casson.C
+++ b/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;
 }
--- a/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/scalarTransport
+++ b/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;
--- a/tutorials/heatTransfer/chtMultiRegionFoam/windshieldCondensation/system/controlDict
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/windshieldCondensation/system/controlDict
@ -61,8 +61,6 @@ functions
        fvOptions
        {
        }
-
-        writeControl    writeTime;
    }

    fileUpdate
--- a/wmake/rules/General/cgal
+++ b/wmake/rules/General/cgal
@ -1,26 +1,21 @@
 # ----------------------------------------------------------------------------
 # CGAL definitions - several possibilities
 #
-#  - missing
-#  - header-only
-#  - header-only, no mpfr
-#  - library, no mpfr
-#  - library, with mpfr (default for older CGAL)
+#  0. missing
+#  1. header-only
+#  2. library, no mpfr
+#  3. library, with mpfr (a likely default)
 #
 # Dispatch according to the defined 'CGAL_FLAVOUR'
 # - names may change [see wmake/scripts/have_cgal]
-# (no-cgal | cgal-header | cgal-header-no-mpfr | cgal-no-mpfr | cgal-mpfr)
+# (no-cgal | cgal-header | 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

 # ----------------------------------------------------------------------------
--- a/wmake/rules/General/cgal-header-no-mpfr
+++ b/wmake/rules/General/cgal-header-no-mpfr
@ -1,18 +0,0 @@
-# -----------------------------------------------------------------------------
-# 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
-
-# -----------------------------------------------------------------------------
--- a/wmake/rules/darwin64Clang/cgal
+++ b/wmake/rules/darwin64Clang/cgal
@ -0,0 +1,14 @@
+# ----------------------------------------------------------------------------
+# 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
+
+# ----------------------------------------------------------------------------
--- a/wmake/rules/linux64Mingw/c++
+++ b/wmake/rules/linux64Mingw/c++
@ -35,7 +35,6 @@ 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
Author	SHA1	Message	Date
alonzameret	44c76c5806	Changed code data structures to correspond with openfoam objects	2023-02-27 12:07:41 +02:00
alonzameret	aaf0245795	Changed code data structures to correspond with openfoam objects	2023-02-27 10:04:47 +02:00
alonzameret	69df0ad468	ENH: new renumbering method based on Hamiltonian path - contributed by Alon Zameret	2023-01-31 16:25:57 +01:00
alonzameret	fdfd2d15a9	Added new renumbering method based on Hamiltonian path	2023-01-17 13:51:52 +02:00