Changed code data structures to correspond with openfoam objects

- contributed by Alon Zameret

META-INFO/api-info

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

doc/Build.md

@@ -1,23 +1,23 @@
 ## Getting the code
 
-Links to all code packs are available on https://dl.openfoam.com. For OpenFOAM-v2206:
+Links to all code packs are available on https://dl.openfoam.com. For OpenFOAM-v2212:
 
 - https://dl.openfoam.com/source/latest/
-- Source: https://dl.openfoam.com/source/v2206/OpenFOAM-v2206.tgz
-- ThirdParty: https://dl.openfoam.com/source/v2206/ThirdParty-v2206.tgz
+- Source: https://dl.openfoam.com/source/v2212/OpenFOAM-v2212.tgz
+- ThirdParty: https://dl.openfoam.com/source/v2212/ThirdParty-v2212.tgz
 
 ## OpenFOAM® Quick Build Guide
 
 Prior to building, ensure that the [system requirements][link openfoam-require]
 are satisfied (including any special [cross-compiling][wiki-cross-compile]
 considerations), and source the correct OpenFOAM environment.
-For example, for the OpenFOAM-v2206 version:
+For example, for the OpenFOAM-v2212 version:
 ```
-source <installation path>/OpenFOAM-v2206/etc/bashrc
+source <installation path>/OpenFOAM-v2212/etc/bashrc
 ```
 e.g. if installed under the `~/openfoam` directory
 ```
-source ~/openfoam/OpenFOAM-v2206/etc/bashrc
+source ~/openfoam/OpenFOAM-v2212/etc/bashrc
 ```
 
 

src/finiteArea/Make/files

@@ -55,7 +55,6 @@ $(basicFaPatchFields)/transform/transformFaPatchScalarField.C
 constraintFaPatchFields = $(faPatchFields)/constraint
 $(constraintFaPatchFields)/empty/emptyFaPatchFields.C
 $(constraintFaPatchFields)/processor/processorFaPatchFields.C
-$(constraintFaPatchFields)/processor/processorFaPatchScalarField.C
 $(constraintFaPatchFields)/wedge/wedgeFaPatchFields.C
 $(constraintFaPatchFields)/wedge/wedgeFaPatchScalarField.C
 $(constraintFaPatchFields)/cyclic/cyclicFaPatchFields.C

src/finiteArea/fields/faPatchFields/basic/coupled/coupledFaPatchField.H

@@ -181,13 +181,6 @@ public:
 
         // Coupled interface functionality
 
-            //- Transform given patch component field
-            virtual void transformCoupleField
-            (
-                solveScalarField& f,
-                const direction cmpt
-            ) const = 0;
-
             //- Update result field based on interface functionality
             virtual void updateInterfaceMatrix
             (

src/finiteArea/fields/faPatchFields/basic/transform/transformFaPatchFieldsFwd.H

@@ -1,58 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2016-2017 Wikki Ltd
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef transformFaPatchFieldsFwd_H
-#define transformFaPatchFieldsFwd_H
-
-#include "fieldTypes.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-template<class Type> class transformFaPatchField;
-
-// typedef transformFaPatchField<scalar> transformFaPatchScalarField;
-// typedef transformFaPatchField<vector> transformFaPatchVectorField;
-// typedef transformFaPatchField<tensor> transformFaPatchTensorField;
-
-// template<class Type> class transformFaPatchField;
-
-makeFaPatchTypeFieldTypedefs(transform)
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/finiteArea/fields/faPatchFields/constraint/cyclic/cyclicFaPatchField.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -187,20 +187,7 @@ void Foam::cyclicFaPatchField<Type>::updateInterfaceMatrix
     transformCoupleField(pnf, cmpt);
 
     // Multiply the field by coefficients and add into the result
-    if (add)
-    {
-        forAll(faceCells, elemI)
-        {
-            result[faceCells[elemI]] += coeffs[elemI]*pnf[elemI];
-        }
-    }
-    else
-    {
-        forAll(faceCells, elemI)
-        {
-            result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
-        }
-    }
+    this->addToInternalField(result, !add, faceCells, coeffs, pnf);
 }
 
 
@@ -227,21 +214,11 @@ void Foam::cyclicFaPatchField<Type>::updateInterfaceMatrix
         pnf[facei + sizeby2] = psiInternal[faceCells[facei]];
     }
 
+    // Transform according to the transformation tensors
+    transformCoupleField(pnf);
+
     // Multiply the field by coefficients and add into the result
-    if (add)
-    {
-        forAll(faceCells, elemI)
-        {
-            result[faceCells[elemI]] += coeffs[elemI]*pnf[elemI];
-        }
-    }
-    else
-    {
-        forAll(faceCells, elemI)
-        {
-            result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
-        }
-    }
+    this->addToInternalField(result, !add, faceCells, coeffs, pnf);
 }
 
 

src/finiteArea/fields/faPatchFields/constraint/cyclic/cyclicFaPatchField.H

@@ -164,16 +164,6 @@ public:
 
         // Coupled interface functionality
 
-            //- Transform neighbour field
-            virtual void transformCoupleField
-            (
-                solveScalarField& f,
-                const direction cmpt
-            ) const
-            {
-                cyclicLduInterfaceField::transformCoupleField(f, cmpt);
-            }
-
             //- Update result field based on interface functionality
             virtual void updateInterfaceMatrix
             (

src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchField.C

@@ -130,13 +130,6 @@ Foam::processorFaPatchField<Type>::processorFaPatchField
 {}
 
 
-// * * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * //
-
-template<class Type>
-Foam::processorFaPatchField<Type>::~processorFaPatchField()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
@@ -203,6 +196,8 @@ void Foam::processorFaPatchField<Type>::initInterfaceMatrixUpdate
         commsType,
         this->patch().patchInternalField(psiInternal)()
     );
+
+    const_cast<processorFaPatchField<Type>&>(*this).updatedMatrix() = false;
 }
 
 
@@ -219,32 +214,28 @@ void Foam::processorFaPatchField<Type>::updateInterfaceMatrix
     const Pstream::commsTypes commsType
 ) const
 {
+    if (this->updatedMatrix())
+    {
+        return;
+    }
+
+    const labelUList& faceCells = this->patch().edgeFaces();
+
     solveScalarField pnf
     (
-        procPatch_.receive<solveScalar>(commsType, this->size())()
+        procPatch_.receive<solveScalar>(commsType, this->size())
     );
 
-    // Transform according to the transformation tensor
-    transformCoupleField(pnf, cmpt);
+    if (!std::is_arithmetic<Type>::value)
+    {
+        // Transform non-scalar data according to the transformation tensor
+        transformCoupleField(pnf, cmpt);
+    }
 
     // Multiply the field by coefficients and add into the result
+    this->addToInternalField(result, !add, faceCells, coeffs, pnf);
 
-    const labelUList& edgeFaces = this->patch().edgeFaces();
-
-    if (add)
-    {
-        forAll(edgeFaces, elemI)
-        {
-            result[edgeFaces[elemI]] += coeffs[elemI]*pnf[elemI];
-        }
-    }
-    else
-    {
-        forAll(edgeFaces, elemI)
-        {
-            result[edgeFaces[elemI]] -= coeffs[elemI]*pnf[elemI];
-        }
-    }
+    const_cast<processorFaPatchField<Type>&>(*this).updatedMatrix() = true;
 }
 
 
@@ -265,6 +256,8 @@ void Foam::processorFaPatchField<Type>::initInterfaceMatrixUpdate
         commsType,
         this->patch().patchInternalField(psiInternal)()
     );
+
+    const_cast<processorFaPatchField<Type>&>(*this).updatedMatrix() = false;
 }
 
 
@@ -280,29 +273,25 @@ void Foam::processorFaPatchField<Type>::updateInterfaceMatrix
     const Pstream::commsTypes commsType
 ) const
 {
+    if (this->updatedMatrix())
+    {
+        return;
+    }
+
+    const labelUList& faceCells = this->patch().edgeFaces();
+
     Field<Type> pnf
     (
-        procPatch_.receive<Type>(commsType, this->size())()
+        procPatch_.receive<Type>(commsType, this->size())
     );
 
+    // Transform according to the transformation tensor
+    transformCoupleField(pnf);
+
     // Multiply the field by coefficients and add into the result
+    this->addToInternalField(result, !add, faceCells, coeffs, pnf);
 
-    const labelUList& edgeFaces = this->patch().edgeFaces();
-
-    if (add)
-    {
-        forAll(edgeFaces, elemI)
-        {
-            result[edgeFaces[elemI]] += coeffs[elemI]*pnf[elemI];
-        }
-    }
-    else
-    {
-        forAll(edgeFaces, elemI)
-        {
-            result[edgeFaces[elemI]] -= coeffs[elemI]*pnf[elemI];
-        }
-    }
+    const_cast<processorFaPatchField<Type>&>(*this).updatedMatrix() = true;
 }
 
 

src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchField.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -38,8 +38,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef processorFaPatchField_H
-#define processorFaPatchField_H
+#ifndef Foam_processorFaPatchField_H
+#define Foam_processorFaPatchField_H
 
 #include "coupledFaPatchField.H"
 #include "processorLduInterfaceField.H"
@@ -139,9 +139,8 @@ public:
         }
 
 
-    // Destructor
-
-        ~processorFaPatchField();
+    //- Destructor
+    ~processorFaPatchField() = default;
 
 
     // Member functions
@@ -169,17 +168,8 @@ public:
             //- Return patch-normal gradient
             virtual tmp<Field<Type>> snGrad() const;
 
-        // Coupled interface functionality
 
-            //- Transform neighbour field
-            virtual void transformCoupleField
-            (
-                solveScalarField& f,
-                const direction cmpt
-            ) const
-            {
-                processorLduInterfaceField::transformCoupleField(f, cmpt);
-            }
+        // Coupled interface functionality
 
             //- Initialise neighbour matrix update
             virtual void initInterfaceMatrixUpdate

src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchFields.H

@@ -25,10 +25,10 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef processorFaPatchFields_H
-#define processorFaPatchFields_H
+#ifndef Foam_processorFaPatchFields_H
+#define Foam_processorFaPatchFields_H
 
-#include "processorFaPatchScalarField.H"
+#include "processorFaPatchField.H"
 #include "fieldTypes.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchFieldsFwd.H

@@ -1,52 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2016-2017 Wikki Ltd
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef processorFaPatchFieldsFwd_H
-#define processorFaPatchFieldsFwd_H
-
-#include "fieldTypes.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-template<class Type> class processorFaPatchField;
-
-makeFaPatchTypeFieldTypedefs(processor)
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchScalarField.C

@@ -1,100 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2019 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "processorFaPatchScalarField.H"
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<>
-void Foam::processorFaPatchField<Foam::scalar>::transformCoupleField
-(
-    solveScalarField& f,
-    const direction cmpt
-) const
-{}
-
-
-template<>
-void Foam::processorFaPatchField<Foam::scalar>::initInterfaceMatrixUpdate
-(
-    solveScalarField& result,
-    const bool add,
-    const lduAddressing& lduAddr,
-    const label patchId,
-    const solveScalarField& psiInternal,
-    const scalarField& coeffs,
-    const direction,
-    const Pstream::commsTypes commsType
-) const
-{
-    procPatch_.send
-    (
-        commsType,
-        patch().patchInternalField(psiInternal)()
-    );
-}
-
-
-template<>
-void Foam::processorFaPatchField<Foam::scalar>::updateInterfaceMatrix
-(
-    solveScalarField& result,
-    const bool add,
-    const lduAddressing& lduAddr,
-    const label patchId,
-    const solveScalarField&,
-    const scalarField& coeffs,
-    const direction,
-    const Pstream::commsTypes commsType
-) const
-{
-    solveScalarField pnf
-    (
-        procPatch_.receive<solveScalar>(commsType, this->size())()
-    );
-
-    const labelUList& edgeFaces = patch().edgeFaces();
-
-    if (add)
-    {
-        forAll(edgeFaces, facei)
-        {
-            result[edgeFaces[facei]] += coeffs[facei]*pnf[facei];
-        }
-    }
-    else
-    {
-        forAll(edgeFaces, facei)
-        {
-            result[edgeFaces[facei]] -= coeffs[facei]*pnf[facei];
-        }
-    }
-}
-
-
-// ************************************************************************* //

src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchScalarField.H

@@ -1,86 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2019 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef processorFaPatchScalarField_H
-#define processorFaPatchScalarField_H
-
-#include "processorFaPatchField.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<>
-void processorFaPatchField<scalar>::transformCoupleField
-(
-    solveScalarField& f,
-    const direction cmpt
-) const;
-
-
-template<>
-void processorFaPatchField<scalar>::initInterfaceMatrixUpdate
-(
-    solveScalarField& result,
-    const bool add,
-    const lduAddressing& lduAddr,
-    const label patchId,
-    const solveScalarField&,
-    const scalarField& coeffs,
-    const direction,
-    const Pstream::commsTypes commsType
-) const;
-
-
-template<>
-void processorFaPatchField<scalar>::updateInterfaceMatrix
-(
-    solveScalarField& result,
-    const bool add,
-    const lduAddressing& lduAddr,
-    const label patchId,
-    const solveScalarField&,
-    const scalarField& coeffs,
-    const direction,
-    const Pstream::commsTypes commsType
-) const;
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/fvOptions/corrections/limitVelocity/limitVelocity.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -71,7 +71,7 @@ Foam::fv::limitVelocity::limitVelocity
 :
     fv::cellSetOption(name, modelType, dict, mesh),
     writeFile(mesh, name, typeName, dict, false),
-    UName_(word::null),
+    UName_("U"),
     max_(0)
 {
     read(dict);

src/meshTools/coordSet/writers/vtk/foamVtkCoordSetWriter.C

@@ -80,6 +80,19 @@ void Foam::vtk::coordSetWriter::beginPiece()
         }
     }
 
+    // Update sizes, similar to
+    // vtk::polyWriter::beginPiece(const pointField&, const edgeList&)
+
+    numberOfPoints_ = nLocalPoints_;
+    numberOfCells_  = nLocalLines_;
+
+    // if (parallel_)
+    // {
+    //     reduce(numberOfPoints_, sumOp<label>());
+    //     reduce(numberOfCells_,  sumOp<label>());
+    // }
+
+
     // Nothing else to do for legacy
     if (legacy()) return;
 
@@ -88,7 +101,7 @@ void Foam::vtk::coordSetWriter::beginPiece()
         format().openTag
         (
             vtk::fileTag::PIECE,
-            vtk::fileAttr::NUMBER_OF_POINTS, nLocalPoints_
+            vtk::fileAttr::NUMBER_OF_POINTS, numberOfPoints_
         );
         if (nLocalVerts_)
         {
@@ -105,7 +118,7 @@ void Foam::vtk::coordSetWriter::beginPiece()
 
 void Foam::vtk::coordSetWriter::writePoints()
 {
-    this->beginPoints(nLocalPoints_);
+    this->beginPoints(numberOfPoints_);  //<- same as nLocalPoints_
 
     {
         for (const pointField& pts : points_)

src/overset/cellCellStencil/cellCellStencil/cellCellStencil.C

@@ -411,7 +411,7 @@ void Foam::cellCellStencil::setUpFrontOnOversetPatch
 void Foam::cellCellStencil::walkFront
 (
     const globalIndex& globalCells,
-    const label layerRelax,
+    const scalar layerRelax,
     const labelListList& allStencil,
     labelList& allCellTypes,
     scalarField& allWeight,
@@ -611,9 +611,11 @@ void Foam::cellCellStencil::walkFront
                                     allWeightWork[nei] = fraction[facei];
                                     allCellTypesWork[nei] = INTERPOLATED;
 
-                                    const label donorId =  compactStencil[nei][0];
+                                    const label donorId =
+                                        compactStencil[nei][0];
 
-                                    volRatio[nei] = V[nei]/compactCellVol[donorId];
+                                    volRatio[nei] =
+                                        V[nei]/compactCellVol[donorId];
 
                                     seedCell
                                     (

src/overset/cellCellStencil/cellCellStencil/cellCellStencil.H

@@ -276,7 +276,7 @@ public:
         void walkFront
         (
             const globalIndex& globalCells,
-            const label layerRelax,
+            const scalar layerRelax,
             const labelListList& allStencil,
             labelList& allCellTypes,
             scalarField& allWeight,

src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.C

@@ -1798,7 +1798,7 @@ Foam::cellCellStencils::inverseDistance::~inverseDistance()
 
 bool Foam::cellCellStencils::inverseDistance::update()
 {
-    label layerRelax(dict_.getOrDefault("layerRelax", 1));
+    scalar layerRelax(dict_.getOrDefault("layerRelax", 1.0));
 
     scalar tol = dict_.getOrDefault("tolerance", 1e-10);
     smallVec_ = mesh_.bounds().span()*tol;

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

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;
+    }
+}
+
+
+// ************************************************************************* //

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

tutorials/incompressible/simpleFoam/rotatingCylinders/0/U

@@ -0,0 +1,39 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2212                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    innerWall 
+    {
+        type            noSlip;
+    }
+
+    outerWall
+    {
+        type            noSlip;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //

tutorials/incompressible/simpleFoam/rotatingCylinders/0/p

@@ -0,0 +1,34 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2212                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    "innerWall|outerWall" 
+    {
+        type            zeroGradient;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //

tutorials/incompressible/simpleFoam/rotatingCylinders/Allrun

@@ -0,0 +1,10 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+runApplication blockMesh
+runApplication $(getApplication)
+./plot
+
+#------------------------------------------------------------------------------

tutorials/incompressible/simpleFoam/rotatingCylinders/constant/MRFProperties

@@ -0,0 +1,30 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2212                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      MRFProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+MRF1
+{
+    cellZone    all;
+    active      yes;
+
+    nonRotatingPatches (outerWall);
+
+    origin    (0 0 0);
+    axis      (0 0 1);
+    omega     100;
+}
+
+// ************************************************************************* //

tutorials/incompressible/simpleFoam/rotatingCylinders/constant/transportProperties

@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2212                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+transportModel  Newtonian;
+
+nu              1;
+
+// ************************************************************************* //

tutorials/incompressible/simpleFoam/rotatingCylinders/constant/turbulenceProperties

@@ -0,0 +1,20 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2212                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType laminar;
+
+// ************************************************************************* //

tutorials/incompressible/simpleFoam/rotatingCylinders/plot

@@ -0,0 +1,42 @@
+#!/bin/sh
+# Require gnuplot
+command -v gnuplot >/dev/null || {
+    echo "FOAM FATAL ERROR: gnuplot not found - skipping graph creation" 1>&2
+    exit 1
+}
+
+
+gnuplot<<EOF
+    set terminal pngcairo font "helvetica,16" size 800,600
+    set xlabel "Radius, r"
+    set xrange [1:2]
+    set grid
+
+    Omega1 = 100.
+    Omega2 = 0.
+    R1 = 1.
+    R2 = 2.
+
+    mu = Omega2/Omega1
+
+    A = Omega1*(1. - R2*R2*mu/R1/R1)/(1. - R2*R2/R1/R1)
+    B = R1*R1*Omega1*(1. - mu)/(1. - R1*R1/R2/R2)
+
+    analytical_utheta(r) = A*r + B/r
+    C = 0.5*A*A*R1*R1 + 2.*A*B*log(R1) - 0.5*B*B/R1/R1
+    analytical_p(r) = 0.5*A*A*r*r + 2.*A*B*log(r) - 0.5*B*B/r/r - C
+
+    set ylabel "Pressure, p"
+    set output "p.png"
+    set samples 20
+    set key bottom right
+    plot "postProcessing/sample1/1000/centreLine_p.xy" w lines lw 2 lc "black" t "OpenFOAM", \
+        analytical_p(x) w p ps 2 pt 6 lw 2 lc rgb "red" t "analytical"
+
+    set ylabel "U_{/Symbol q}, rad/s"
+    set output "Utheta.png"
+    set samples 20
+    set key top right
+    plot "postProcessing/sample1/1000/centreLine_U_U:Transformed.xy" u 1:6 w lines lw 2 lc "black" t "OpenFOAM", \
+        analytical_utheta(x) w p ps 2 pt 6 lw 2 lc rgb "red" t "analytical"
+EOF

tutorials/incompressible/simpleFoam/rotatingCylinders/system/blockMeshDict

@@ -0,0 +1,126 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2212                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale 1;
+
+geom
+{
+    r0          1;
+    r1          2;
+
+    nr          40;
+    ntheta      40;
+
+    g           10;
+    
+    mr0         #eval{ - $r0 };
+    mr1         #eval{ - $r1 };
+    invG        #eval{ 1./$g };
+}
+
+vertices
+(
+    ( $:geom.r0         0 0)
+    (         0  $geom.r0 0)
+    ($:geom.mr0         0 0)
+    (         0 $geom.mr0 0)
+    ( $:geom.r1         0 0)
+    (         0  $geom.r1 0)
+    ($:geom.mr1         0 0)
+    (         0 $geom.mr1 0)
+
+    ( $:geom.r0         0 1)
+    (         0  $geom.r0 1)
+    ($:geom.mr0         0 1)
+    (         0 $geom.mr0 1)
+    ( $:geom.r1         0 1)
+    (         0  $geom.r1 1)
+    ($:geom.mr1         0 1)
+    (         0 $geom.mr1 1)
+);
+
+blockInfo
+all
+($:geom.ntheta $:geom.nr 1)
+simpleGrading (1 ((0.5 0.5 $:geom.g)(0.5 0.5 $:geom.invG)) 1);
+
+blocks
+(
+    hex (1 0 4 5  9  8 12 13) $blockInfo
+    hex (2 1 5 6 10  9 13 14) $blockInfo
+    hex (3 2 6 7 11 10 14 15) $blockInfo
+    hex (0 3 7 4  8 11 15 12) $blockInfo
+);
+
+edges
+(
+    arc 0 1 origin (0 0 0)
+    arc 1 2 origin (0 0 0)
+    arc 2 3 origin (0 0 0)
+    arc 3 0 origin (0 0 0)
+    arc 8 9 origin (0 0 1)
+    arc 9 10 origin (0 0 1)
+    arc 10 11 origin (0 0 1)
+    arc 11 8 origin (0 0 1)
+
+    arc 4 5 origin (0 0 0)
+    arc 5 6 origin (0 0 0)
+    arc 6 7 origin (0 0 0)
+    arc 7 4 origin (0 0 0)
+    arc 12 13 origin (0 0 1)
+    arc 13 14 origin (0 0 1)
+    arc 14 15 origin (0 0 1)
+    arc 15 12 origin (0 0 1)
+);
+
+defaultPatch
+{
+    name    frontAndBack;
+    type    empty;
+}
+
+boundary
+(
+    innerWall
+    {
+        type wall;
+        faces
+        (
+            (1 0  8  9)
+            (2 1  9 10)
+            (3 2 10 11)
+            (0 3 11 8)
+        );
+    }
+
+    outerWall
+    {
+        type wall;
+        faces
+        (
+            (5 13 12 4)
+            (6 14 13 5)
+            (7 15 14 6)
+            (4 12 15 7)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //

tutorials/incompressible/simpleFoam/rotatingCylinders/system/controlDict

@@ -0,0 +1,93 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2212                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     simpleFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         1000;
+
+deltaT          1;
+
+writeControl    timeStep;
+
+writeInterval   10;
+
+purgeWrite      3;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable true;
+
+functions
+{
+    transform1
+    {
+        type            fieldCoordinateSystemTransform;
+        libs            (fieldFunctionObjects);
+        writeControl    writeTime;
+
+        fields          (U);
+        coordinateSystem
+        {
+            origin          (0 0 0);
+            coordinateRotation
+            {
+                type            axesRotation;
+                e3              (0 0 1);
+                e1              (1 0 0);
+            }
+        }
+    }
+
+    sample1
+    {
+        type            sets;
+        libs            (sampling);
+        writeControl    writeTime;
+        fields          (U p U:Transformed);
+        interpolationScheme cellPoint;
+        setFormat       raw;
+
+        sets
+        (
+            centreLine
+            {
+                type        uniform;
+                axis        x;
+                start       (1 0 0);
+                end         (2 0 0);
+                nPoints     20;
+            }
+        );
+    }
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/simpleFoam/rotatingCylinders/system/fvSchemes

@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2212                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         steadyState;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+    div(phi,U)      bounded Gauss linearUpwind grad(U);
+    div((nuEff*dev2(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/simpleFoam/rotatingCylinders/system/fvSolution

@@ -0,0 +1,57 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2212                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p
+    {
+        solver          GAMG;
+        tolerance       1e-10;
+        relTol          0.1;
+        smoother        GaussSeidel;
+    }
+
+    U
+    {
+        solver          smoothSolver;
+        smoother        GaussSeidel;
+        tolerance       1e-10;
+        relTol          0.1;
+    }
+}
+
+SIMPLE
+{
+    nNonOrthogonalCorrectors 0;
+    pRefCell        0;
+    pRefValue       0;
+}
+
+relaxationFactors
+{
+    fields
+    {
+        p               0.3;
+    }
+    equations
+    {
+        U               0.7;
+    }
+}
+
+
+// ************************************************************************* //