From 9891d371725e42a981af55eec39ec84dd433ea9a Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Wed, 8 Oct 2008 07:55:07 +0100
Subject: [PATCH] parallel postChannel

---
 .../miscellaneous/postChannel/Make/options    |   2 +
 .../miscellaneous/postChannel/channelIndex.C  | 381 ++++++++++++------
 .../miscellaneous/postChannel/channelIndex.H  | 114 ++++--
 .../postChannel/channelIndexTemplates.C       | 101 +++++
 .../miscellaneous/postChannel/collapse.H      |   8 +-
 .../miscellaneous/postChannel/postChannel.C   |  15 +-
 .../miscellaneous/postChannel/postChannelDict |  14 +-
 .../miscellaneous/postChannel/sumData.C       |  49 +++
 .../miscellaneous/postChannel/sumData.H       | 200 +++++++++
 .../miscellaneous/postChannel/sumDataI.H      | 227 +++++++++++
 .../channel395/constant/postChannelDict       |  14 +-
 11 files changed, 930 insertions(+), 195 deletions(-)
 create mode 100644 applications/utilities/postProcessing/miscellaneous/postChannel/channelIndexTemplates.C
 create mode 100644 applications/utilities/postProcessing/miscellaneous/postChannel/sumData.C
 create mode 100644 applications/utilities/postProcessing/miscellaneous/postChannel/sumData.H
 create mode 100644 applications/utilities/postProcessing/miscellaneous/postChannel/sumDataI.H

diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/Make/options b/applications/utilities/postProcessing/miscellaneous/postChannel/Make/options
index d38cd8b180..b90b6fdd4e 100644
--- a/applications/utilities/postProcessing/miscellaneous/postChannel/Make/options
+++ b/applications/utilities/postProcessing/miscellaneous/postChannel/Make/options
@@ -1,7 +1,9 @@
 EXE_INC = \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude
 
 EXE_LIBS = \
+    -lmeshTools \
     -lfiniteVolume \
     -lsampling
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.C b/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.C
index 74315de31a..cdba8c0201 100644
--- a/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.C
+++ b/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.C
@@ -25,156 +25,269 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "channelIndex.H"
+#include "boolList.H"
+#include "syncTools.H"
+#include "OFstream.H"
+#include "meshTools.H"
+#include "Time.H"
+#include "SortableList.H"
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * * * //
 
-channelIndex::channelIndex(const fvMesh& m)
-:
-    indexingDict_
-    (
-        IOobject
-        (
-            "postChannelDict",
-            m.time().constant(),
-            m,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    ),
-    nx_(readLabel(indexingDict_.lookup("Nx"))),
-    ny_(indexingDict_.lookup("Ny")),
-    nz_(readLabel(indexingDict_.lookup("Nz"))),
-    symmetric_
-    (
-        readBool(indexingDict_.lookup("symmetric"))
-    ),
-    cumNy_(ny_.size()),
-    nLayers_(ny_[0])
+template<>
+const char* Foam::NamedEnum<Foam::vector::components, 3>::names[] =
 {
-    // initialise the layers
-    cumNy_[0] = ny_[0];
+    "x",
+    "y",
+    "z"
+};
 
-    for (label j=1; j<ny_.size(); j++)
+const Foam::NamedEnum<Foam::vector::components, 3>
+    Foam::channelIndex::vectorComponentsNames_;
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Determines face blocking
+void Foam::channelIndex::walkOppositeFaces
+(
+    const polyMesh& mesh,
+    const labelList& startFaces,
+    boolList& blockedFace
+)
+{
+    const cellList& cells = mesh.cells();
+    const faceList& faces = mesh.faces();
+    label nBnd = mesh.nFaces() - mesh.nInternalFaces();
+
+    DynamicList<label> frontFaces(startFaces);
+    forAll(frontFaces, i)
     {
-        nLayers_ += ny_[j];
-        cumNy_[j] = ny_[j]+cumNy_[j-1];
+        label faceI = frontFaces[i];
+        blockedFace[faceI] = true;
+    }
+
+    while (returnReduce(frontFaces.size(), sumOp<label>()) > 0)
+    {
+        // Transfer across.
+        boolList isFrontBndFace(nBnd, false);
+        forAll(frontFaces, i)
+        {
+            label faceI = frontFaces[i];
+
+            if (!mesh.isInternalFace(faceI))
+            {
+                isFrontBndFace[faceI-mesh.nInternalFaces()] = true;
+            }
+        }
+        syncTools::swapBoundaryFaceList(mesh, isFrontBndFace, false);
+
+        // Add 
+        forAll(isFrontBndFace, i)
+        {
+            label faceI = mesh.nInternalFaces()+i;
+            if (isFrontBndFace[i] && !blockedFace[faceI])
+            {
+                blockedFace[faceI] = true;
+                frontFaces.append(faceI);
+            }
+        }
+
+        // Transfer across cells
+        DynamicList<label> newFrontFaces(frontFaces.size());
+
+        forAll(frontFaces, i)
+        {
+            label faceI = frontFaces[i];
+
+            {
+                const cell& ownCell = cells[mesh.faceOwner()[faceI]];
+
+                label oppositeFaceI = ownCell.opposingFaceLabel(faceI, faces);
+
+                if (oppositeFaceI == -1)
+                {
+                    FatalErrorIn("channelIndex::walkOppositeFaces(..)")
+                        << "Face:" << faceI << " owner cell:" << ownCell
+                        << " is not a hex?" << abort(FatalError);
+                }
+                else
+                {
+                    if (!blockedFace[oppositeFaceI])
+                    {
+                        blockedFace[oppositeFaceI] = true;
+                        newFrontFaces.append(oppositeFaceI);
+                    }
+                }
+            }
+
+            if (mesh.isInternalFace(faceI))
+            {
+                const cell& neiCell = mesh.cells()[mesh.faceNeighbour()[faceI]];
+
+                label oppositeFaceI = neiCell.opposingFaceLabel(faceI, faces);
+
+                if (oppositeFaceI == -1)
+                {
+                    FatalErrorIn("channelIndex::walkOppositeFaces(..)")
+                        << "Face:" << faceI << " neighbour cell:" << neiCell
+                        << " is not a hex?" << abort(FatalError);
+                }
+                else
+                {
+                    if (!blockedFace[oppositeFaceI])
+                    {
+                        blockedFace[oppositeFaceI] = true;
+                        newFrontFaces.append(oppositeFaceI);
+                    }
+                }
+            }
+        }
+
+        frontFaces.transfer(newFrontFaces);
     }
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+// Calculate regions.
+void Foam::channelIndex::calcLayeredRegions
+(
+    const polyMesh& mesh,
+    const labelList& startFaces
+)
+{
+    boolList blockedFace(mesh.nFaces(), false);
+    walkOppositeFaces
+    (
+        mesh,
+        startFaces,
+        blockedFace
+    );
 
-channelIndex::~channelIndex()
-{}
+
+    if (false)
+    {
+        OFstream str(mesh.time().path()/"blockedFaces.obj");
+        label vertI = 0;
+        forAll(blockedFace, faceI)
+        {
+            if (blockedFace[faceI])
+            {
+                const face& f = mesh.faces()[faceI];
+                forAll(f, fp)
+                {
+                    meshTools::writeOBJ(str, mesh.points()[f[fp]]);
+                }
+                str<< 'f';
+                forAll(f, fp)
+                {
+                    str << ' ' << vertI+fp+1;
+                }
+                str << nl;
+                vertI += f.size();
+            }
+        }
+    }
+
+
+    // Do analysis for connected regions
+    cellRegion_.reset(new regionSplit(mesh, blockedFace));
+
+    Info<< "Detected " << cellRegion_().nRegions() << " layers." << nl << endl;
+
+    // Sum number of entries per region
+    regionCount_ = regionSum(scalarField(mesh.nCells(), 1.0));
+
+    // Average cell centres to determine ordering.
+    pointField regionCc
+    (
+        regionSum(mesh.cellCentres())
+      / regionCount_
+    );
+
+    SortableList<scalar> sortComponent(regionCc.component(dir_));
+
+    sortMap_ = sortComponent.indices();
+
+    y_ = sortComponent;
+
+    if (symmetric_)
+    {
+        y_.setSize(cellRegion_().nRegions()/2);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::channelIndex::channelIndex
+(
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    symmetric_(readBool(dict.lookup("symmetric"))),
+    dir_(vectorComponentsNames_.read(dict.lookup("component")))
+{
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    const wordList patchNames(dict.lookup("patches"));
+
+    label nFaces = 0;
+
+    forAll(patchNames, i)
+    {
+        label patchI = patches.findPatchID(patchNames[i]);
+
+        if (patchI == -1)
+        {
+            FatalErrorIn("channelIndex::channelIndex(const polyMesh&)")
+                << "Illegal patch " << patchNames[i]
+                << ". Valid patches are " << patches.name()
+                << exit(FatalError);
+        }
+
+        nFaces += patches[patchI].size();
+    }
+
+    labelList startFaces(nFaces);
+    nFaces = 0;
+
+    forAll(patchNames, i)
+    {
+        const polyPatch& pp = patches[patches.findPatchID(patchNames[i])];
+
+        forAll(pp, j)
+        {
+            startFaces[nFaces++] = pp.start()+j;
+        }
+    }
+
+    // Calculate regions.
+    calcLayeredRegions(mesh, startFaces);
+}
+
+
+Foam::channelIndex::channelIndex
+(
+    const polyMesh& mesh,
+    const labelList& startFaces,
+    const bool symmetric,
+    const direction dir
+)
+:
+    symmetric_(symmetric),
+    dir_(dir)
+{
+    // Calculate regions.
+    calcLayeredRegions(mesh, startFaces);
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-scalarField channelIndex::collapse
-(
-    const volScalarField& vsf,
-    const bool asymmetric
-) const
-{
-    scalarField cs(nLayers(), 0.0);
-
-    forAll(cs, j)
-    {
-        // sweep over all cells in this layer
-        for (label i=0; i<nx(); i++)
-        {
-            for (label k=0; k<nz(); k++)
-            {
-                cs[j] += vsf[operator()(i,j,k)];
-            }
-        }
-
-        // and divide by the number of cells in the layer
-        cs[j] /= scalar(nx()*nz());
-    }
-
-    if (symmetric_)
-    {
-        label nlb2 = nLayers()/2;
-
-        if (asymmetric)
-        {
-            for (label j=0; j<nlb2; j++)
-            {
-                cs[j] = 0.5*(cs[j] - cs[nLayers() - j - 1]);
-            }
-        }
-        else
-        {
-            for (label j=0; j<nlb2; j++)
-            {
-                cs[j] = 0.5*(cs[j] + cs[nLayers() - j - 1]);
-            }
-        }
-
-        cs.setSize(nlb2);
-    }
-
-    return cs;
-}
-
-
-scalarField channelIndex::y
-(
-    const volVectorField& cellCentres
-) const
-{
-    if (symmetric_)
-    {
-        scalarField Y(nLayers()/2);
-
-        for (label j=0; j<nLayers()/2; j++)
-        {
-            Y[j] = cellCentres[operator()(0, j, 0)].y();
-        }
-
-        return Y;
-    }
-    else
-    {
-        scalarField Y(nLayers());
-
-        for (label j=0; j<nLayers(); j++)
-        {
-            Y[j] = cellCentres[operator()(0, j, 0)].y();
-        }
-
-        return Y;
-    }
-}
-
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
-label channelIndex::operator()
-(
-    const label Jx,
-    const label Jy,
-    const label Jz
-) const
-{
-    label index(0);
-
-    // count up `full' layers in the mesh
-    label j(0);
-    label tmpJy(Jy);
-
-    while(Jy >= cumNy_[j])
-    {
-        index += nx_*ny_[j]*nz_;
-        tmpJy -= ny_[j];
-        j++;
-    }
-
-    index += Jx + nx_*tmpJy + nx_*ny_[j]*Jz;
-
-    return index;
-}
-
 
 // ************************************************************************* //
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.H b/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.H
index 500d13f4eb..796e20160c 100644
--- a/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.H
+++ b/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.H
@@ -26,19 +26,25 @@ Class
     Foam::channelIndex
 
 Description
-    does indexing for a standard (Christer-style) channel.
-Assumes that the blocks are aranged in the y direction.
+    Does averaging of fields over layers of cells. Assumes layered mesh.
 
 SourceFiles
     channelIndex.C
-    channelIndexIO.C
 
 \*---------------------------------------------------------------------------*/
 
 #ifndef channelIndex_H
 #define channelIndex_H
 
-#include "fvCFD.H"
+#include "regionSplit.H"
+#include "direction.H"
+#include "scalarField.H"
+#include "polyMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
 
 
 /*---------------------------------------------------------------------------*\
@@ -47,21 +53,47 @@ SourceFiles
 
 class channelIndex
 {
+
     // Private data
 
-        IOdictionary indexingDict_;
+        static const NamedEnum<vector::components, 3> vectorComponentsNames_;
 
-        const label nx_;
-        const labelList ny_;
-        const label nz_;
+        //- Is mesh symmetric
         const bool symmetric_;
 
-        labelList cumNy_;
-        label nLayers_;
+        //- direction to sort
+        const direction dir_;
+
+        //- Per cell the global region
+        autoPtr<regionSplit> cellRegion_;
+
+        //- Per global region the number of cells (scalarField so we can use
+        //  field algebra)
+        scalarField regionCount_;
+
+        //- From sorted region back to unsorted global region
+        labelList sortMap_;
+
+        //- Sorted component of cell centres
+        scalarField y_;
+
 
 
     // Private Member Functions
 
+        void walkOppositeFaces
+        (
+            const polyMesh& mesh,
+            const labelList& startFaces,
+            boolList& blockedFace
+        );
+
+        void calcLayeredRegions
+        (
+            const polyMesh& mesh,
+            const labelList& startFaces
+        );
+
         //- Disallow default bitwise copy construct and assignment
         channelIndex(const channelIndex&);
         void operator=(const channelIndex&);
@@ -71,56 +103,54 @@ public:
 
     // Constructors
 
-        channelIndex(const fvMesh& m);
+        //- Construct from dictionary
+        channelIndex(const polyMesh&, const dictionary&);
 
-
-    // Destructor
-
-        ~channelIndex();
+        //- Construct from supplied starting faces
+        channelIndex
+        (
+            const polyMesh& mesh,
+            const labelList& startFaces,
+            const bool symmetric,
+            const direction dir
+        );
 
 
     // Member Functions
 
         // Access
 
-            //- return number of layers
-            label nLayers() const
-            {
-                return nLayers_;
-            }
-
-            //- number of cells in X direction
-            label nx() const
-            {
-                return nx_;
-            }
-
-            //- number of cells in Z direction
-            label nz() const
-            {
-                return nz_;
-            }
+            //- Sum field per region
+            template<class T>
+            Field<T> regionSum(const Field<T>& cellField) const;
 
             //- collapse a field to a line
-            scalarField collapse
+            template<class T>
+            Field<T> collapse
             (
-                const volScalarField& vsf,
+                const Field<T>& vsf,
                 const bool asymmetric=false
             ) const;
 
             //- return the field of Y locations from the cell centres
-            scalarField y
-            (
-                const volVectorField& cellCentres
-            ) const;
+            const scalarField& y() const
+            {
+                return y_;
+            }
 
-
-    // Member Operators
-
-        label operator()(const label, const label, const label) const;
 };
 
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "channelIndexTemplates.C"
+#endif
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndexTemplates.C b/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndexTemplates.C
new file mode 100644
index 0000000000..9fc8ff0360
--- /dev/null
+++ b/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndexTemplates.C
@@ -0,0 +1,101 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "channelIndex.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class T>
+Foam::Field<T> Foam::channelIndex::regionSum(const Field<T>& cellField) const
+{
+    Field<T> regionField(cellRegion_().nRegions(), pTraits<T>::zero);
+
+    forAll(cellRegion_(), cellI)
+    {
+        regionField[cellRegion_()[cellI]] += cellField[cellI];
+    }
+
+    // Global sum
+    Pstream::listCombineGather(regionField, plusEqOp<T>());
+    Pstream::listCombineScatter(regionField);
+
+    return regionField;
+}
+
+
+template<class T>
+Foam::Field<T> Foam::channelIndex::collapse
+(
+    const Field<T>& cellField,
+    const bool asymmetric
+) const
+{
+    // Average and order
+    Field<T> regionField
+    (
+        regionSum(cellField)
+      / regionCount_,
+        sortMap_
+    );
+
+    // Symmetry?
+    if (symmetric_)
+    {
+        label nlb2 = cellRegion_().nRegions()/2;
+
+        if (asymmetric)
+        {
+            for (label j=0; j<nlb2; j++)
+            {
+                regionField[j] =
+                    0.5
+                  * (
+                        regionField[j]
+                      - regionField[cellRegion_().nRegions() - j - 1]
+                    );
+            }
+        }
+        else
+        {
+            for (label j=0; j<nlb2; j++)
+            {
+                regionField[j] =
+                    0.5
+                  * (
+                        regionField[j]
+                      + regionField[cellRegion_().nRegions() - j - 1]
+                    );
+            }
+        }
+
+        regionField.setSize(nlb2);
+    }
+
+    return regionField;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/collapse.H b/applications/utilities/postProcessing/miscellaneous/postChannel/collapse.H
index ec047e5862..b3bf559411 100644
--- a/applications/utilities/postProcessing/miscellaneous/postChannel/collapse.H
+++ b/applications/utilities/postProcessing/miscellaneous/postChannel/collapse.H
@@ -1,16 +1,16 @@
     scalarField UMeanXvalues = channelIndexing.collapse
     (
-        UMean.component(vector::X)
+        UMean.component(vector::X)()
     );
 
     scalarField UMeanYvalues = channelIndexing.collapse
     (
-        UMean.component(vector::Y)
+        UMean.component(vector::Y)()
     );
 
     scalarField UMeanZvalues = channelIndexing.collapse
     (
-        UMean.component(vector::Z)
+        UMean.component(vector::Z)()
     );
 
     scalarField RxxValues = channelIndexing.collapse(Rxx);
@@ -38,7 +38,7 @@
         0.5*(sqr(urmsValues) + sqr(vrmsValues) + sqr(wrmsValues));
 
 
-    scalarField y = channelIndexing.y(mesh.C());
+    const scalarField& y = channelIndexing.y();
 
     makeGraph(y, UMeanXvalues, "Uf", UMean.path(), gFormat);
     makeGraph(y, urmsValues, "u", UMean.path(), gFormat);
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/postChannel.C b/applications/utilities/postProcessing/miscellaneous/postChannel/postChannel.C
index 2de4c534a8..08051cb1da 100644
--- a/applications/utilities/postProcessing/miscellaneous/postChannel/postChannel.C
+++ b/applications/utilities/postProcessing/miscellaneous/postChannel/postChannel.C
@@ -67,7 +67,20 @@ int main(int argc, char *argv[])
     const word& gFormat = runTime.graphFormat();
 
     // Setup channel indexing for averaging over channel down to a line
-    channelIndex channelIndexing(mesh);
+
+    IOdictionary channelDict
+    (
+        IOobject
+        (
+            "postChannelDict",
+            mesh.time().constant(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+    channelIndex channelIndexing(mesh, channelDict);
+
 
     // For each time step read all fields
     for (label i=startTime; i<endTime; i++)
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/postChannelDict b/applications/utilities/postProcessing/miscellaneous/postChannel/postChannelDict
index e69c0386ef..d63cd65667 100644
--- a/applications/utilities/postProcessing/miscellaneous/postChannel/postChannelDict
+++ b/applications/utilities/postProcessing/miscellaneous/postChannel/postChannelDict
@@ -15,14 +15,14 @@ FoamFile
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-Nx 40; 
-Ny 
-( 
- 25 
- 25 
-); 
-Nz 30; 
+// Seed patches to start layering from
+patches (bottomWall);
 
+// Direction in which the layers are
+component y;
+
+// Is the mesh symmetric? If so average(symmetric fields) or
+// subtract(asymmetric) contributions from both halves
 symmetric true;
 
 // ************************************************************************* //
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/sumData.C b/applications/utilities/postProcessing/miscellaneous/postChannel/sumData.C
new file mode 100644
index 0000000000..8662ec578b
--- /dev/null
+++ b/applications/utilities/postProcessing/miscellaneous/postChannel/sumData.C
@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "sumData.H"
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+    Foam::Ostream& os,
+    const Foam::sumData& wDist
+)
+{
+    return os
+        << wDist.oldFace_ << token::SPACE
+        << wDist.sum_ << token::SPACE << wDist.count_;
+}
+
+
+Foam::Istream& Foam::operator>>(Foam::Istream& is, Foam::sumData& wDist)
+{
+    return is >> wDist.oldFace_ >> wDist.sum_ >> wDist.count_;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/sumData.H b/applications/utilities/postProcessing/miscellaneous/postChannel/sumData.H
new file mode 100644
index 0000000000..15a26a3aaa
--- /dev/null
+++ b/applications/utilities/postProcessing/miscellaneous/postChannel/sumData.H
@@ -0,0 +1,200 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::sumData
+
+Description
+    Sums data while walking across cells. Used in collapsing fields.
+
+SourceFiles
+    sumDataI.H
+    sumData.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef sumData_H
+#define sumData_H
+
+#include "point.H"
+#include "tensor.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class polyPatch;
+class polyMesh;
+
+/*---------------------------------------------------------------------------*\
+                           Class sumData Declaration
+\*---------------------------------------------------------------------------*/
+
+class sumData
+{
+    // Private data
+
+        //- Previous face
+        label oldFace_;
+
+        //- summed data
+        scalar sum_;
+
+        //- number of items summed
+        label count_;
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        inline sumData();
+
+        //- Construct from count
+        inline sumData
+        (
+            const label oldFace,
+            const scalar sum,
+            const label count
+        );
+
+
+    // Member Functions
+
+        // Access
+
+            inline label oldFace() const
+            {
+                return oldFace_;
+            }
+
+            inline scalar sum() const
+            {
+                return sum_;
+            }
+
+            inline label count() const
+            {
+                return count_;
+            }
+
+
+        // Needed by FaceCellWave
+
+            //- Check whether origin has been changed at all or
+            //  still contains original (invalid) value.
+            inline bool valid() const;
+
+            //- Check for identical geometrical data. Used for cyclics checking.
+            inline bool sameGeometry
+            (
+                const polyMesh&,
+                const sumData&,
+                const scalar
+            ) const;
+
+            //- Convert any absolute coordinates into relative to (patch)face
+            //  centre
+            inline void leaveDomain
+            (
+                const polyMesh&,
+                const polyPatch&,
+                const label patchFaceI,
+                const point& faceCentre
+            );
+
+            //- Reverse of leaveDomain
+            inline void enterDomain
+            (
+                const polyMesh&,
+                const polyPatch&,
+                const label patchFaceI,
+                const point& faceCentre
+            );
+
+            //- Apply rotation matrix to any coordinates
+            inline void transform
+            (
+                const polyMesh&,
+                const tensor&
+            );
+
+            //- Influence of neighbouring face.
+            inline bool updateCell
+            (
+                const polyMesh&,
+                const label thisCellI,
+                const label neighbourFaceI,
+                const sumData& neighbourInfo,
+                const scalar tol
+            );
+
+            //- Influence of neighbouring cell.
+            inline bool updateFace
+            (
+                const polyMesh&,
+                const label thisFaceI,
+                const label neighbourCellI,
+                const sumData& neighbourInfo,
+                const scalar tol
+            );
+
+            //- Influence of different value on same face.
+            inline bool updateFace
+            (
+                const polyMesh&,
+                const label thisFaceI,
+                const sumData& neighbourInfo,
+                const scalar tol
+            );
+
+    // Member Operators
+
+        // Needed for List IO
+        inline bool operator==(const sumData&) const;
+
+        inline bool operator!=(const sumData&) const;
+
+
+    // IOstream Operators
+
+        friend Ostream& operator<<(Ostream&, const sumData&);
+        friend Istream& operator>>(Istream&, sumData&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "sumDataI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/sumDataI.H b/applications/utilities/postProcessing/miscellaneous/postChannel/sumDataI.H
new file mode 100644
index 0000000000..107cba063d
--- /dev/null
+++ b/applications/utilities/postProcessing/miscellaneous/postChannel/sumDataI.H
@@ -0,0 +1,227 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "polyMesh.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Null constructor
+inline Foam::sumData::sumData()
+:
+    oldFace_(-1),
+    sum_(0.0),
+    count_(0)
+{}
+
+
+// Construct from components
+inline Foam::sumData::sumData
+(
+    const label oldFace,
+    const scalar sum,
+    const label count
+)
+:
+    oldFace_(oldFace),
+    sum_(sum),
+    count_(count)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline bool Foam::sumData::valid() const
+{
+    return oldFace_ != -1;
+}
+
+
+// No geometric data so never any problem on cyclics
+inline bool Foam::sumData::sameGeometry
+(
+    const polyMesh&,
+    const sumData&,
+    const scalar
+) const
+{
+    return true;
+}
+
+
+// No geometric data.
+inline void Foam::sumData::leaveDomain
+(
+    const polyMesh&,
+    const polyPatch& patch,
+    const label patchFaceI,
+    const point& faceCentre
+)
+{}
+
+
+// No geometric data.
+inline void Foam::sumData::transform
+(
+    const polyMesh&,
+    const tensor& rotTensor
+)
+{}
+
+
+// No geometric data.
+inline void Foam::sumData::enterDomain
+(
+    const polyMesh&,
+    const polyPatch& patch,
+    const label patchFaceI,
+    const point& faceCentre
+)
+{
+    oldFace_ = -1;
+}
+
+
+// Update cell with neighbouring face information
+inline bool Foam::sumData::updateCell
+(
+    const polyMesh&,
+    const label thisCellI,
+    const label neighbourFaceI,
+    const sumData& neighbourInfo,
+    const scalar tol
+)
+{
+    if (!valid())
+    {
+        FatalErrorIn("sumData::updateCell") << "problem"
+            << abort(FatalError);
+        return false;
+    }
+
+
+    if (count_ == 0)
+    {
+        sum_ += neighbourInfo.sum();
+        count_ = neighbourInfo.count() + 1;
+        oldFace_ = neighbourFaceI;
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// Update face with neighbouring cell information
+inline bool Foam::sumData::updateFace
+(
+    const polyMesh& mesh,
+    const label thisFaceI,
+    const label neighbourCellI,
+    const sumData& neighbourInfo,
+    const scalar tol
+)
+{
+    // From cell to its faces.
+
+    // Check that face is opposite the previous one.
+    const cell& cFaces = mesh.cells()[neighbourCellI];
+
+    label wantedFaceI = cFaces.opposingFaceLabel
+    (
+        neighbourInfo.oldFace(),
+        mesh.faces()
+    );
+
+    if (thisFaceI == wantedFaceI)
+    {
+        if (count_ != 0)
+        {
+            FatalErrorIn("sumData::updateFace") << "problem"
+                << abort(FatalError);
+            return false;
+        }
+
+        sum_ += neighbourInfo.sum();
+        count_ = neighbourInfo.count();
+        oldFace_ = thisFaceI;
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// Update face with coupled face information
+inline bool Foam::sumData::updateFace
+(
+    const polyMesh&,
+    const label thisFaceI,
+    const sumData& neighbourInfo,
+    const scalar tol
+)
+{
+    // From face to face (e.g. coupled faces)
+    if (count_ == 0)
+    {
+        sum_ += neighbourInfo.sum();
+        count_ = neighbourInfo.count();
+        oldFace_ = thisFaceI;
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}    
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+inline bool Foam::sumData::operator==(const Foam::sumData& rhs)
+ const
+{
+    return
+        oldFace() == rhs.oldFace()
+     && sum() == rhs.sum()
+     && count() == rhs.count();
+}
+
+
+inline bool Foam::sumData::operator!=(const Foam::sumData& rhs)
+ const
+{
+    return !(*this == rhs);
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/channelOodles/channel395/constant/postChannelDict b/tutorials/channelOodles/channel395/constant/postChannelDict
index e69c0386ef..d63cd65667 100644
--- a/tutorials/channelOodles/channel395/constant/postChannelDict
+++ b/tutorials/channelOodles/channel395/constant/postChannelDict
@@ -15,14 +15,14 @@ FoamFile
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-Nx 40; 
-Ny 
-( 
- 25 
- 25 
-); 
-Nz 30; 
+// Seed patches to start layering from
+patches (bottomWall);
 
+// Direction in which the layers are
+component y;
+
+// Is the mesh symmetric? If so average(symmetric fields) or
+// subtract(asymmetric) contributions from both halves
 symmetric true;
 
 // ************************************************************************* //