From d5820be68b1874f4ea77631875bf2bb828068efa Mon Sep 17 00:00:00 2001
From: Will Bainbridge <http://cfd.direct>
Date: Thu, 21 Apr 2022 12:25:25 +0100
Subject: [PATCH] layerAverage: Re-write layer generation using mesh wave

This change means this function is determining the sequence in which
points are plotted topologically. This makes it possible to plot a layer
average along a pipe that goes through many changes of direction.

Previously, the function determined the order by means of a geometric
sort in the plot direction. This only worked when the layers were
perpendicular to one of the coordinate axes.
---
 .../postProcessing/graphs/graphLayerAverage   |   6 +-
 .../meshes/meshShapes/cell/oppositeCellFace.C |  21 +-
 .../field/layerAverage/layerAverage.C         | 281 ++++++------------
 .../field/layerAverage/layerAverage.H         |  25 +-
 .../layerAverage/layerAverageTemplates.C      |   6 +-
 .../field/layerAverage/layerInfo.H            | 200 +++++++++++++
 .../field/layerAverage/layerInfoI.H           | 244 +++++++++++++++
 7 files changed, 559 insertions(+), 224 deletions(-)
 create mode 100644 src/functionObjects/field/layerAverage/layerInfo.H
 create mode 100644 src/functionObjects/field/layerAverage/layerInfoI.H

diff --git a/etc/caseDicts/postProcessing/graphs/graphLayerAverage b/etc/caseDicts/postProcessing/graphs/graphLayerAverage
index c114db9a79..23223ac58a 100644
--- a/etc/caseDicts/postProcessing/graphs/graphLayerAverage
+++ b/etc/caseDicts/postProcessing/graphs/graphLayerAverage
@@ -14,8 +14,10 @@ patches         (<patchNames>);  // Patches from which layers extrude
 
 zones           ();              // Zones from which layers extrude
 
-axis            <axisName>;      // The independent variable of the graph. Can
-                                 // be "x", "y" or "z".
+axis            distance;        // The independent variable of the graph. Can
+                                 // be "x", "y", "z", "xyz" (all coordinates
+                                 // written out), or "distance" (from the start
+                                 // point).
 
 symmetric       false;           // Are the layers symmetric about the centre?
 
diff --git a/src/OpenFOAM/meshes/meshShapes/cell/oppositeCellFace.C b/src/OpenFOAM/meshes/meshShapes/cell/oppositeCellFace.C
index 24f6a6a63b..72b42de90f 100644
--- a/src/OpenFOAM/meshes/meshShapes/cell/oppositeCellFace.C
+++ b/src/OpenFOAM/meshes/meshShapes/cell/oppositeCellFace.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -42,11 +42,13 @@ Foam::label Foam::cell::opposingFaceLabel
     const faceUList& meshFaces
 ) const
 {
-    // Algorithm:
     // Go through all the faces of the cell and find the one which
-    // does not share a single vertex with the master face.  If there
-    // are two or more such faces, return the first one and issue a
-    // warning; if there is no opposite face, return -1;
+    // does not share a single vertex with the master face
+    //
+    // If there are no such faces (e.g., a tetrahedron), return -1
+    //
+    // If there are more than one such faces (e.g., a hex with split faces),
+    // return -2
 
     const face& masterFace = meshFaces[masterFaceLabel];
 
@@ -96,12 +98,9 @@ Foam::label Foam::cell::opposingFaceLabel
                 }
                 else
                 {
-                    // There has already been an opposite face.
-                    // Non-prismatic cell
-                    Info<< "Multiple faces not sharing vertex: "
-                        << oppositeFaceLabel << " and "
-                        << curFaceLabels[facei] << endl;
-                    return -1;
+                    // There has already been a face with no shared points.
+                    // This cell is not prismatic.
+                    return -2;
                 }
             }
         }
diff --git a/src/functionObjects/field/layerAverage/layerAverage.C b/src/functionObjects/field/layerAverage/layerAverage.C
index 51a98df1d9..e322ef4d87 100644
--- a/src/functionObjects/field/layerAverage/layerAverage.C
+++ b/src/functionObjects/field/layerAverage/layerAverage.C
@@ -24,6 +24,8 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "layerAverage.H"
+#include "FaceCellWave.H"
+#include "layerInfo.H"
 #include "regionSplit.H"
 #include "syncTools.H"
 #include "volFields.H"
@@ -39,204 +41,61 @@ namespace Foam
         defineTypeNameAndDebug(layerAverage, 0);
         addToRunTimeSelectionTable(functionObject, layerAverage, dictionary);
     }
-
-    template<>
-    const char*
-        Foam::NamedEnum<Foam::vector::components, 3>::names[] =
-        {"x", "y", "z"};
 }
 
-const Foam::NamedEnum<Foam::vector::components, 3>
-    Foam::functionObjects::layerAverage::axisNames_;
-
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::functionObjects::layerAverage::walkOppositeFaces
-(
-    const labelList& startFaces,
-    const boolList& startFaceIntoOwners,
-    boolList& blockedFace,
-    boolList& cellIsLayer
-)
-{
-    // Initialise the front
-    DynamicList<label> frontFaces(startFaces);
-    DynamicList<bool> frontFaceIntoOwners(startFaceIntoOwners);
-    DynamicList<label> newFrontFaces(frontFaces.size());
-    DynamicList<bool> newFrontFaceIntoOwners(frontFaceIntoOwners.size());
-
-    // Set the start and end faces as blocked
-    UIndirectList<bool>(blockedFace, startFaces) = true;
-
-    // Iterate until the front is empty
-    while (returnReduce(frontFaces.size(), sumOp<label>()) > 0)
-    {
-        // Transfer front faces across couplings
-        boolList bndFaceIsFront(mesh_.nFaces() - mesh_.nInternalFaces(), false);
-        forAll(frontFaces, i)
-        {
-            const label facei = frontFaces[i];
-            const label intoOwner = frontFaceIntoOwners[i];
-
-            if (!mesh_.isInternalFace(facei) && !intoOwner)
-            {
-                bndFaceIsFront[facei - mesh_.nInternalFaces()] = true;
-            }
-        }
-        syncTools::swapBoundaryFaceList(mesh_, bndFaceIsFront);
-
-        // Set new front faces
-        forAll(bndFaceIsFront, bndFacei)
-        {
-            const label facei = mesh_.nInternalFaces() + bndFacei;
-
-            if (bndFaceIsFront[bndFacei] && !blockedFace[facei])
-            {
-                blockedFace[facei] = true;
-                frontFaces.append(facei);
-                frontFaceIntoOwners.append(true);
-            }
-        }
-
-        // Transfer across cells
-        newFrontFaces.clear();
-        newFrontFaceIntoOwners.clear();
-
-        forAll(frontFaces, i)
-        {
-            const label facei = frontFaces[i];
-            const label intoOwner = frontFaceIntoOwners[i];
-
-            const label celli =
-                intoOwner
-              ? mesh_.faceOwner()[facei]
-              : mesh_.isInternalFace(facei) ? mesh_.faceNeighbour()[facei] : -1;
-
-            if (celli != -1)
-            {
-                cellIsLayer[celli] = true;
-
-                const label oppositeFacei =
-                    mesh_.cells()[celli]
-                   .opposingFaceLabel(facei, mesh_.faces());
-                const bool oppositeIntoOwner =
-                    mesh_.faceOwner()[oppositeFacei] != celli;
-
-                if (oppositeFacei == -1)
-                {
-                    FatalErrorInFunction
-                        << "Cannot find face on cell " << mesh_.cells()[celli]
-                        << " opposing face " << mesh_.faces()[facei]
-                        << ". Mesh is not layered." << exit(FatalError);
-                }
-                else if (!blockedFace[oppositeFacei])
-                {
-                    blockedFace[oppositeFacei] = true;
-                    newFrontFaces.append(oppositeFacei);
-                    newFrontFaceIntoOwners.append(oppositeIntoOwner);
-                }
-            }
-        }
-
-        frontFaces.transfer(newFrontFaces);
-        frontFaceIntoOwners.transfer(newFrontFaceIntoOwners);
-    }
-
-    // Determine whether cells on the other sides of couplings are in layers
-    boolList bndCellIsLayer(mesh_.nFaces() - mesh_.nInternalFaces(), false);
-    forAll(bndCellIsLayer, bndFacei)
-    {
-        const label facei = mesh_.nInternalFaces() + bndFacei;
-        const label owni = mesh_.faceOwner()[facei];
-
-        bndCellIsLayer[bndFacei] = cellIsLayer[owni];
-    }
-    syncTools::swapBoundaryFaceList(mesh_, bndCellIsLayer);
-
-    // Block faces where one adjacent cell is in a layer and the other is not
-    forAll(mesh_.faces(), facei)
-    {
-        const label owni = mesh_.faceOwner()[facei];
-        if (mesh_.isInternalFace(facei))
-        {
-            const label nbri = mesh_.faceNeighbour()[facei];
-            if (cellIsLayer[owni] != cellIsLayer[nbri])
-            {
-                blockedFace[facei] = true;
-            }
-        }
-        else
-        {
-            const label bndFacei = facei - mesh_.nInternalFaces();
-            if (cellIsLayer[owni] != bndCellIsLayer[bndFacei])
-            {
-                blockedFace[facei] = true;
-            }
-        }
-    }
-}
-
-
 void Foam::functionObjects::layerAverage::calcLayers()
 {
+    // Initialise the faces from which the layers extrude
     DynamicList<label> startFaces;
-    DynamicList<bool> startFaceIntoOwners;
-
+    DynamicList<layerInfo> startFacesInfo;
     forAll(patchIDs_, i)
     {
         const polyPatch& pp = mesh_.boundaryMesh()[patchIDs_[i]];
-        startFaces.append(identity(pp.size()) + pp.start());
-        startFaceIntoOwners.append(boolList(pp.size(), true));
+        forAll(pp, j)
+        {
+            startFaces.append(pp.start() + j);
+            startFacesInfo.append(layerInfo(0, -1));
+        }
     }
-
     forAll(zoneIDs_, i)
     {
         const faceZone& fz = mesh_.faceZones()[zoneIDs_[i]];
-        startFaces.append(fz);
-        startFaceIntoOwners.append(fz.flipMap());
+        forAll(fz, j)
+        {
+            startFaces.append(fz[j]);
+            startFacesInfo.append(layerInfo(0, fz.flipMap()[j] ? -1 : +1));
+        }
     }
 
-    // Identify faces that separate the layers
-    boolList blockedFace(mesh_.nFaces(), false);
-    boolList cellIsLayer(mesh_.nCells(), false);
-    walkOppositeFaces
+    // Wave to generate layer indices
+    List<layerInfo> allFaceLayerInfo(mesh_.nFaces());
+    List<layerInfo> allCellLayerInfo(mesh_.nCells());
+    FaceCellWave<layerInfo> wave
     (
+        mesh_,
         startFaces,
-        startFaceIntoOwners,
-        blockedFace,
-        cellIsLayer
+        startFacesInfo,
+        allFaceLayerInfo,
+        allCellLayerInfo,
+        mesh_.globalData().nTotalCells() + 1
     );
 
-    // Do analysis for connected layers
-    regionSplit rs(mesh_, blockedFace);
-    nLayers_ = rs.nRegions();
-    cellLayer_.transfer(rs);
-
-    // Get rid of regions that are not layers
-    label layeri0 = labelMax, layeri1 = -labelMax;
+    // Copy indices out of the wave and determine the total number of layers
+    nLayers_ = 0;
+    cellLayer_ = labelList(mesh_.nCells(), -1);
     forAll(cellLayer_, celli)
     {
-        if (cellIsLayer[celli])
+        if (allCellLayerInfo[celli].valid(wave.data()))
         {
-            layeri0 = min(layeri0, cellLayer_[celli]);
-            layeri1 = max(layeri1, cellLayer_[celli]);
-        }
-        else
-        {
-            cellLayer_[celli] = -1;
-        }
-    }
-    reduce(layeri0, minOp<label>());
-    reduce(layeri1, maxOp<label>());
-    nLayers_ = layeri0 != labelMax ? 1 + layeri1 - layeri0 : 0;
-    forAll(cellLayer_, celli)
-    {
-        if (cellLayer_[celli] != -1)
-        {
-            cellLayer_[celli] -= layeri0;
+            const label layeri = allCellLayerInfo[celli].cellLayer();
+            nLayers_ = max(nLayers_, layeri + 1);
+            cellLayer_[celli] = layeri;
         }
     }
+    reduce(nLayers_, maxOp<label>());
 
     // Report
     if (nLayers_ != 0)
@@ -245,26 +104,33 @@ void Foam::functionObjects::layerAverage::calcLayers()
     }
     else
     {
-        WarningInFunction
-            << "No layers detected" << endl;
+        WarningInFunction<< "No layers detected" << endl;
+    }
+
+    // Write the indices for debugging
+    if (debug)
+    {
+        tmp<volScalarField> cellLayer =
+            volScalarField::New
+            (
+                "cellLayer",
+                mesh_,
+                dimensionedScalar(dimless, 0)
+            );
+        cellLayer.ref().primitiveFieldRef() = List<scalar>(cellLayer_);
+        cellLayer.ref().write();
     }
 
     // Sum number of entries per layer
     layerCount_ = sum(scalarField(mesh_.nCells(), 1));
 
     // Average the cell centres
-    const pointField layerCentres(sum(mesh_.cellCentres())/layerCount_);
-
-    // Sort by the direction component
-    x_ = layerCentres.component(axis_);
-    sortMap_ = identity(layerCentres.size());
-    stableSort(sortMap_, UList<scalar>::less(x_));
-    x_.map(x_, sortMap_);
+    layerCentre_ = sum(mesh_.cellCentres())/layerCount_;
 
     // If symmetric, keep only half of the coordinates
     if (symmetric_)
     {
-        x_.setSize(nLayers_/2);
+        layerCentre_.setSize(nLayers_/2);
     }
 }
 
@@ -273,8 +139,28 @@ template<>
 Foam::vector
 Foam::functionObjects::layerAverage::symmetricCoeff<Foam::vector>() const
 {
+    direction i = -1;
+
+    switch (axis_)
+    {
+        case coordSet::axisType::X:
+        case coordSet::axisType::Y:
+        case coordSet::axisType::Z:
+            i = label(axis_) - label(coordSet::axisType::X);
+            break;
+        case coordSet::axisType::XYZ:
+        case coordSet::axisType::DISTANCE:
+        case coordSet::axisType::DEFAULT:
+            FatalErrorInFunction
+                << "Symmetric layer average requested with "
+                << coordSet::axisTypeNames_[axis_] << " axis. Symmetric "
+                << "averaging is only possible along coordinate axes."
+                << exit(FatalError);
+            break;
+    }
+
     vector result = vector::one;
-    result[axis_] = -1;
+    result[i] = -1;
     return result;
 }
 
@@ -343,7 +229,15 @@ bool Foam::functionObjects::layerAverage::read(const dictionary& dict)
 
     symmetric_ = dict.lookupOrDefault<bool>("symmetric", false);
 
-    axis_ = axisNames_.read(dict.lookup("axis"));
+    axis_ =
+        coordSet::axisTypeNames_
+        [
+            dict.lookupOrDefault<word>
+            (
+                "axis",
+                coordSet::axisTypeNames_[coordSet::axisType::DEFAULT]
+            )
+        ];
 
     fields_ = dict.lookup<wordList>("fields");
 
@@ -422,13 +316,28 @@ bool Foam::functionObjects::layerAverage::write()
     }
 
     // Write
-    if (Pstream::master() && x_.size())
+    if (Pstream::master() && layerCentre_.size())
     {
+        scalarField layerDistance(layerCentre_.size(), 0);
+        for (label i = 1; i < layerCentre_.size(); ++ i)
+        {
+            layerDistance[i] =
+                layerDistance[i-1] + mag(layerCentre_[i] - layerCentre_[i-1]);
+        }
+
         formatter_->write
         (
             outputPath/mesh_.time().timeName(),
             typeName,
-            coordSet(true, axisNames_[axis_], x_),
+            coordSet
+            (
+                identity(layerCentre_.size()),
+                word::null,
+                layerCentre_,
+                coordSet::axisTypeNames_[coordSet::axisType::DISTANCE],
+                layerDistance,
+                coordSet::axisTypeNames_[axis_]
+            ),
             fieldNames
             #define TypeValueSetsParameter(Type, nullArg) , Type##ValueSets
             FOR_ALL_FIELD_TYPES(TypeValueSetsParameter)
diff --git a/src/functionObjects/field/layerAverage/layerAverage.H b/src/functionObjects/field/layerAverage/layerAverage.H
index eff19e227d..489eaa5b5e 100644
--- a/src/functionObjects/field/layerAverage/layerAverage.H
+++ b/src/functionObjects/field/layerAverage/layerAverage.H
@@ -90,12 +90,6 @@ class layerAverage
 :
     public fvMeshFunctionObject
 {
-    // Private Static Data
-
-        //- Names of vector components
-        static const NamedEnum<vector::components, 3> axisNames_;
-
-
     // Private Data
 
         //- Patches which form the start of the layers
@@ -111,7 +105,7 @@ class layerAverage
         bool symmetric_;
 
         //- The direction over which to plot the results
-        vector::components axis_;
+        coordSet::axisType axis_;
 
         //- Per cell the global layer
         label nLayers_;
@@ -122,11 +116,8 @@ class layerAverage
         //- Per global layer the number of cells
         scalarField layerCount_;
 
-        //- From sorted layer back to unsorted global layer
-        labelList sortMap_;
-
-        //- Sorted component of cell centres
-        scalarField x_;
+        //- The average centre of each layer
+        pointField layerCentre_;
 
         //- Fields to sample
         wordList fields_;
@@ -137,16 +128,6 @@ class layerAverage
 
     // Private Member Functions
 
-        //- Walk through layers marking faces that separate layers and cells
-        //  that are within layers
-        void walkOppositeFaces
-        (
-            const labelList& startFaces,
-            const boolList& startFaceIntoOwners,
-            boolList& blockedFace,
-            boolList& cellIsLayer
-        );
-
         //- Create the layer information, the sort map, and the scalar axis
         void calcLayers();
 
diff --git a/src/functionObjects/field/layerAverage/layerAverageTemplates.C b/src/functionObjects/field/layerAverage/layerAverageTemplates.C
index 1c897c5977..27154e3b62 100644
--- a/src/functionObjects/field/layerAverage/layerAverageTemplates.C
+++ b/src/functionObjects/field/layerAverage/layerAverageTemplates.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -63,8 +63,8 @@ Foam::Field<T> Foam::functionObjects::layerAverage::average
     const Field<T>& cellField
 ) const
 {
-    // Sum, average and sort
-    Field<T> layerField(sum(cellField)/layerCount_, sortMap_);
+    // Sum and average
+    Field<T> layerField(sum(cellField)/layerCount_);
 
     // Handle symmetry
     if (symmetric_)
diff --git a/src/functionObjects/field/layerAverage/layerInfo.H b/src/functionObjects/field/layerAverage/layerInfo.H
new file mode 100644
index 0000000000..aa62b0aea5
--- /dev/null
+++ b/src/functionObjects/field/layerAverage/layerInfo.H
@@ -0,0 +1,200 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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::layerInfo
+
+Description
+    Class to be used with FaceCellWave which enumerates layers of cells
+
+SourceFiles
+    layerInfoI.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef layerInfo_H
+#define layerInfo_H
+
+#include "pointField.H"
+#include "face.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class polyPatch;
+class polyMesh;
+class transformer;
+
+// Forward declaration of friend functions and operators
+class layerInfo;
+Ostream& operator<<(Ostream&, const layerInfo&);
+Istream& operator>>(Istream&, layerInfo&);
+
+/*---------------------------------------------------------------------------*\
+                          Class layerInfo Declaration
+\*---------------------------------------------------------------------------*/
+
+class layerInfo
+{
+    // Private Data
+
+        //- Layer index
+        label layer_;
+
+        //- Direction of propagation (+1 or -1) relative to the face normal.
+        //  Only valid when on a face. Takes a value of 0 when in a cell.
+        label direction_;
+
+        //- The previous face index. Only valid when in a cell. Takes a value
+        //  of -1 if on a face.
+        label prevFace_;
+
+
+    // Private Member Functions
+
+        //- Generate an error in the event that opposing waves collided
+        void collide() const;
+
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        inline layerInfo();
+
+        //- Construct given the face layer index and direction
+        inline layerInfo(const label faceLayer, const label direction);
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return the face layer index
+            inline label faceLayer() const;
+
+            //- Return the face layer index
+            inline label cellLayer() const;
+
+
+        // Needed by FaceCellWave
+
+            //- Check whether the layerInfo has been changed at all or still
+            //  contains original (invalid) value.
+            template<class TrackingData>
+            inline bool valid(TrackingData& td) const;
+
+            //- Check for identical geometrical data. Used for checking
+            //  consistency across cyclics.
+            template<class TrackingData>
+            inline bool sameGeometry
+            (
+                const polyMesh&,
+                const layerInfo&,
+                const scalar,
+                TrackingData& td
+            ) const;
+
+            //- Transform across an interface
+            template<class TrackingData>
+            inline void transform
+            (
+                const polyPatch& patch,
+                const label patchFacei,
+                const transformer& transform,
+                TrackingData& td
+            );
+
+            //- Influence of neighbouring face
+            template<class TrackingData>
+            inline bool updateCell
+            (
+                const polyMesh&,
+                const label thisCelli,
+                const label neighbourFacei,
+                const layerInfo& neighbourInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Influence of neighbouring cell
+            template<class TrackingData>
+            inline bool updateFace
+            (
+                const polyMesh&,
+                const label thisFacei,
+                const label neighbourCelli,
+                const layerInfo& neighbourInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Influence of different value on same face
+            template<class TrackingData>
+            inline bool updateFace
+            (
+                const polyMesh&,
+                const label thisFacei,
+                const layerInfo& neighbourInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Test equality
+            template<class TrackingData>
+            inline bool equal
+            (
+                const layerInfo&,
+                TrackingData& td
+            ) const;
+
+
+    // Member Operators
+
+        inline bool operator==(const layerInfo&) const;
+        inline bool operator!=(const layerInfo&) const;
+
+
+    // IOstream Operators
+
+        friend Ostream& operator<<(Ostream&, const layerInfo&);
+        friend Istream& operator>>(Istream&, layerInfo&);
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "layerInfoI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/layerAverage/layerInfoI.H b/src/functionObjects/field/layerAverage/layerInfoI.H
new file mode 100644
index 0000000000..841a17d3b1
--- /dev/null
+++ b/src/functionObjects/field/layerAverage/layerInfoI.H
@@ -0,0 +1,244 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 "layerInfo.H"
+#include "polyMesh.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+inline void Foam::layerInfo::collide() const
+{
+    FatalErrorInFunction
+        << "Layer extrusions collided. Check the patches/zones from which "
+        << "layers are being extruded and ensure that they do not point "
+        << "in opposite directions."
+        << exit(FatalError);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+inline Foam::layerInfo::layerInfo()
+:
+    layer_(-labelMax),
+    direction_(0),
+    prevFace_(-labelMax)
+{}
+
+
+inline Foam::layerInfo::layerInfo(const label faceLayer, const label direction)
+:
+    layer_(2*faceLayer - 1),
+    direction_(direction),
+    prevFace_(-labelMax)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline Foam::label Foam::layerInfo::faceLayer() const
+{
+    if (layer_ % 2 != 1)
+    {
+        FatalError
+            << "Face layer index requested from cell layer info"
+            << exit(FatalError);
+    }
+
+    return (layer_ + 1)/2;
+}
+
+
+inline Foam::label Foam::layerInfo::cellLayer() const
+{
+    if (layer_ % 2 != 0)
+    {
+        FatalError
+            << "Cell layer index requested from face layer info"
+            << exit(FatalError);
+    }
+
+    return layer_/2;
+}
+
+
+template<class TrackingData>
+inline bool Foam::layerInfo::valid(TrackingData& td) const
+{
+    return layer_ != -labelMax;
+}
+
+
+template<class TrackingData>
+inline bool Foam::layerInfo::sameGeometry
+(
+    const polyMesh&,
+    const layerInfo& l,
+    const scalar tol,
+    TrackingData& td
+) const
+{
+    return layer_ == l.layer_;
+}
+
+
+template<class TrackingData>
+inline void Foam::layerInfo::transform
+(
+    const polyPatch& patch,
+    const label patchFacei,
+    const transformer& transform,
+    TrackingData& td
+)
+{}
+
+
+template<class TrackingData>
+inline bool Foam::layerInfo::updateCell
+(
+    const polyMesh& mesh,
+    const label thisCelli,
+    const label neighbourFacei,
+    const layerInfo& neighbourLayerInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    const bool o = thisCelli == mesh.faceOwner()[neighbourFacei];
+
+    if (o == (neighbourLayerInfo.direction_ < 0))
+    {
+        if (valid(td) && prevFace_ != neighbourFacei) collide();
+
+        layer_ = valid(td) ? -labelMax : neighbourLayerInfo.layer_ + 1;
+        direction_ = 0;
+        prevFace_ = neighbourFacei;
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class TrackingData>
+inline bool Foam::layerInfo::updateFace
+(
+    const polyMesh& mesh,
+    const label thisFacei,
+    const label neighbourCelli,
+    const layerInfo& neighbourLayerInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    const cell& c = mesh.cells()[neighbourCelli];
+    const label prevFacei = neighbourLayerInfo.prevFace_;
+    const label nextFacei = c.opposingFaceLabel(prevFacei, mesh.faces());
+
+    if (nextFacei == thisFacei)
+    {
+        const label direction =
+            mesh.faceOwner()[thisFacei] == neighbourCelli ? +1 : -1;
+
+        if (valid(td) && direction != direction_) collide();
+
+        layer_ = valid(td) ? -labelMax : neighbourLayerInfo.layer_ + 1;
+        direction_ = direction;
+        prevFace_ = -labelMax;
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class TrackingData>
+inline bool Foam::layerInfo::updateFace
+(
+    const polyMesh& mesh,
+    const label thisFacei,
+    const layerInfo& neighbourLayerInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    const label direction = -neighbourLayerInfo.direction_;
+
+    if (valid(td) && direction != direction_) collide();
+
+    layer_ = valid(td) ? -labelMax : neighbourLayerInfo.layer_;
+    direction_ = direction;
+    prevFace_ = -labelMax;
+
+    return true;
+}
+
+
+template<class TrackingData>
+inline bool Foam::layerInfo::equal
+(
+    const layerInfo& rhs,
+    TrackingData& td
+) const
+{
+    return operator==(rhs);
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+inline bool Foam::layerInfo::operator==(const Foam::layerInfo& rhs) const
+{
+    return layer_ == rhs.layer_;
+}
+
+
+inline bool Foam::layerInfo::operator!=(const Foam::layerInfo& rhs) const
+{
+    return !(*this == rhs);
+}
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const layerInfo& l)
+{
+    return os << l.layer_ << token::SPACE << l.direction_;
+}
+
+
+Foam::Istream& Foam::operator>>(Istream& is, layerInfo& l)
+{
+    return is >> l.layer_ >> l.direction_;
+}
+
+
+// ************************************************************************* //