From ea2b087db9b8144cdaa5808e650a076f289fc2e3 Mon Sep 17 00:00:00 2001
From: henry <Henry Weller h.weller@opencfd.co.uk>
Date: Tue, 7 Oct 2008 08:28:55 +0100
Subject: [PATCH 01/10] Upgraded to the new improved SVD weighting.

---
 .../schemes/quadraticFit/quadraticFitData.C   | 36 +++++++++++++------
 1 file changed, 25 insertions(+), 11 deletions(-)

diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C
index 0caebcd421..7271bad0ff 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C
@@ -206,7 +206,7 @@ Foam::label Foam::quadraticFitData::calcFit
     // calculate the matrix of the polynomial components
     scalarRectangularMatrix B(C.size(), minSize_, scalar(0));
 
-    for(label ip = 0; ip < C.size(); ip++)
+        for(label ip = 0; ip < C.size(); ip++)
     {
         const point& p = C[ip];
 
@@ -228,9 +228,8 @@ Foam::label Foam::quadraticFitData::calcFit
         pz /= scale;
 
         label is = 0;
-        B[ip][is++] = wts[ip];
-
-        B[ip][is++] = wts[ip]*px;
+        B[ip][is++] = wts[0]*wts[ip];
+        B[ip][is++] = wts[0]*wts[ip]*px;
         B[ip][is++] = wts[ip]*sqr(px);
 
         if (dim_ >= 2)
@@ -255,14 +254,14 @@ Foam::label Foam::quadraticFitData::calcFit
     label nSVDzeros = 0;
 
     bool goodFit = false;
-    for(scalar tol = SMALL; tol < 0.1 && !goodFit; tol *= 10)
+    for(int iIt = 0; iIt < 35 && !goodFit; iIt++)
     {
-        SVD svd(B, tol);
+        SVD svd(B, SMALL);
 
-        scalar fit0 = wts[0]*svd.VSinvUt()[0][0];
-        scalar fit1 = wts[1]*svd.VSinvUt()[0][1];
+        scalar fit0 = wts[0]*wts[0]*svd.VSinvUt()[0][0];
+        scalar fit1 = wts[0]*wts[1]*svd.VSinvUt()[0][1];
 
-        goodFit = sign(fit0) == sign(fit1);
+        goodFit = sign(fit0) == sign(fit1) && fit0 < 1 && fit1 < 1;
 
         if (goodFit)
         {
@@ -270,12 +269,27 @@ Foam::label Foam::quadraticFitData::calcFit
             fit_[faci][1] = fit1;
             for(label i = 2; i < stencilSize; i++)
             {
-                fit_[faci][i] = wts[i]*svd.VSinvUt()[0][i];
+                fit_[faci][i] = wts[0]*wts[i]*svd.VSinvUt()[0][i];
             }
 
             singVals = svd.S();
             nSVDzeros = svd.nZeros();
         }
+        else // (not good fit so increase weight in the centre and for linear)
+        {
+            wts[0] *= 10;
+            wts[1] *= 10;
+            for(label i = 0; i < B.n(); i++)
+            {
+                B[i][0] *= 10;
+                B[i][1] *= 10;
+            }
+            for(label j = 0; j < B.m(); j++)
+            {
+                B[0][j] *= 10;
+                B[1][j] *= 10;
+            }
+        }
     }
     if (!goodFit)
     {
@@ -283,7 +297,7 @@ Foam::label Foam::quadraticFitData::calcFit
         (
             "quadraticFitData::calcFit(const pointField&, const label)"
             ) << "For face " << faci << endl
-            << "Fit not good even once tol >= 0.1"
+            << "Cannot find good fit"
             << exit(FatalError);
     }
 

From 81f08283f4d9bf15e5347ba0889d8f5bbeba62b2 Mon Sep 17 00:00:00 2001
From: henry <Henry Weller h.weller@opencfd.co.uk>
Date: Tue, 7 Oct 2008 10:07:43 +0100
Subject: [PATCH 02/10] Added momentum equation relaxation to allow forcing
 diagonal-dominance.

---
 applications/solvers/incompressible/oodles/oodles.C | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/applications/solvers/incompressible/oodles/oodles.C b/applications/solvers/incompressible/oodles/oodles.C
index eb5fa65dcf..8c51065db1 100644
--- a/applications/solvers/incompressible/oodles/oodles.C
+++ b/applications/solvers/incompressible/oodles/oodles.C
@@ -66,6 +66,9 @@ int main(int argc, char *argv[])
           + sgsModel->divDevBeff(U)
         );
 
+        // Optionally ensure diagonal-dominance of the momentum matrix
+        UEqn.relax();
+
         if (momentumPredictor)
         {
             solve(UEqn == -fvc::grad(p));

From 3c69749b47d6d2a9ecf08bb21e8034df678f86c6 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Tue, 7 Oct 2008 13:10:16 +0100
Subject: [PATCH 03/10] intel compiler from Ayad

---
 .../fieldAverageItem/fieldAverageItem.H       | 32 +++++++++++++------
 1 file changed, 23 insertions(+), 9 deletions(-)

diff --git a/src/postProcessing/fieldAverage/fieldAverageItem/fieldAverageItem.H b/src/postProcessing/fieldAverage/fieldAverageItem/fieldAverageItem.H
index 1bcb5689bb..a8b8830016 100644
--- a/src/postProcessing/fieldAverage/fieldAverageItem/fieldAverageItem.H
+++ b/src/postProcessing/fieldAverage/fieldAverageItem/fieldAverageItem.H
@@ -102,15 +102,6 @@ private:
         baseType base_;
 
 
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-//        fieldAverageItem(const fieldAverageItem&);
-
-        //- Disallow default bitwise assignment
-//        void operator=(const fieldAverageItem&);
-
-
 public:
 
     // Constructors
@@ -175,6 +166,29 @@ public:
 
         void operator=(const fieldAverageItem&);
 
+    // Friend Operators
+
+        friend bool operator==
+        (
+            const fieldAverageItem& a,
+            const fieldAverageItem& b
+        )
+        {
+            return
+                a.fieldName_ == b.fieldName_
+             && a.mean_ == b.mean_
+             && a.prime2Mean_ == b.prime2Mean_
+             && a.base_ == b.base_;
+        }
+
+        friend bool operator!=
+        (
+            const fieldAverageItem& a,
+            const fieldAverageItem& b
+        )
+        {
+            return !(a == b);
+        }
 
     // IOstream Operators
 

From 69c567167d218a7b6e5358cf30bd5b6be9f71de0 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Tue, 7 Oct 2008 18:02:50 +0100
Subject: [PATCH 04/10] intel compiler messages

---
 .../rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.C    | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.C b/applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.C
index 4db4faac50..5a17820680 100644
--- a/applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.C
+++ b/applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.C
@@ -157,9 +157,8 @@ void maxwellSlipUFvPatchVectorField::updateCoeffs()
 
     if(thermalCreep_)
     {
-        const GeometricField<scalar, fvPatchField, volMesh>& vsfT =
-            this->db().objectRegistry::
-            lookupObject<GeometricField<scalar, fvPatchField, volMesh> >("T");
+        const volScalarField& vsfT =
+            this->db().objectRegistry::lookupObject<volScalarField>("T");
         label patchi = this->patch().index();
         const fvPatchScalarField& pT = vsfT.boundaryField()[patchi];
         Field<vector> gradpT = fvc::grad(vsfT)().boundaryField()[patchi];

From 0bbc8fd2c42cf6f8abd27671fe6aff431f493d68 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Tue, 7 Oct 2008 18:24:56 +0100
Subject: [PATCH 05/10] parallel checkMesh

---
 .../manipulation/checkMesh/checkGeometry.C    | 82 +++++++++++------
 .../manipulation/checkMesh/checkTopology.C    | 25 +++--
 .../manipulation/checkMesh/printMeshStats.C   | 91 +++++++++++++------
 .../polyBoundaryMesh/polyBoundaryMesh.C       |  4 +-
 4 files changed, 141 insertions(+), 61 deletions(-)

diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
index 6b1a10d90e..8eb1f3f05f 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
@@ -39,10 +39,15 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
         if (mesh.checkEdgeAlignment(true, validDirs, &nonAlignedPoints))
         {
             noFailedChecks++;
+            label nNonAligned = returnReduce
+            (
+                nonAlignedPoints.size(),
+                sumOp<label>()
+            );
 
-            if (nonAlignedPoints.size() > 0)
+            if (nNonAligned > 0)
             {
-                Pout<< "  <<Writing " << nonAlignedPoints.size()
+                Info<< "  <<Writing " << nNonAligned
                     << " points on non-aligned edges to set "
                     << nonAlignedPoints.name() << endl;
                 nonAlignedPoints.write();
@@ -59,16 +64,21 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
         {
             noFailedChecks++;
 
-            if (cells.size() > 0)
+            label nNonClosed = returnReduce(cells.size(), sumOp<label>());
+
+            if (nNonClosed > 0)
             {
-                Pout<< "  <<Writing " << cells.size()
+                Info<< "  <<Writing " << nNonClosed
                     << " non closed cells to set " << cells.name() << endl;
                 cells.write();
             }
         }
-        if (aspectCells.size() > 0)
+
+        label nHighAspect = returnReduce(aspectCells.size(), sumOp<label>());
+
+        if (nHighAspect > 0)
         {
-            Pout<< "  <<Writing " << aspectCells.size()
+            Info<< "  <<Writing " << nHighAspect
                 << " cells with high aspect ratio to set "
                 << aspectCells.name() << endl;
             aspectCells.write();
@@ -81,9 +91,11 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
         {
             noFailedChecks++;
 
-            if (faces.size() > 0)
+            label nFaces = returnReduce(faces.size(), sumOp<label>());
+
+            if (nFaces > 0)
             {
-                Pout<< "  <<Writing " << faces.size()
+                Info<< "  <<Writing " << nFaces
                     << " zero area faces to set " << faces.name() << endl;
                 faces.write();
             }
@@ -96,9 +108,11 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
         {
             noFailedChecks++;
 
-            if (cells.size() > 0)
+            label nCells = returnReduce(cells.size(), sumOp<label>());
+
+            if (nCells > 0)
             {
-                Pout<< "  <<Writing " << cells.size()
+                Info<< "  <<Writing " << nCells
                     << " zero volume cells to set " << cells.name() << endl;
                 cells.write();
             }
@@ -112,9 +126,11 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
             noFailedChecks++;
         }
 
-        if (faces.size() > 0)
+        label nFaces = returnReduce(faces.size(), sumOp<label>());
+
+        if (nFaces > 0)
         {
-            Pout<< "  <<Writing " << faces.size()
+            Info<< "  <<Writing " << nFaces
                 << " non-orthogonal faces to set " << faces.name() << endl;
             faces.write();
         }
@@ -127,9 +143,11 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
         {
             noFailedChecks++;
 
-            if (faces.size() > 0)
+            label nFaces = returnReduce(faces.size(), sumOp<label>());
+
+            if (nFaces > 0)
             {
-                Pout<< "  <<Writing " << faces.size()
+                Info<< "  <<Writing " << nFaces
                     << " faces with incorrect orientation to set "
                     << faces.name() << endl;
                 faces.write();
@@ -143,9 +161,11 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
         {
             noFailedChecks++;
 
-            if (faces.size() > 0)
+            label nFaces = returnReduce(faces.size(), sumOp<label>());
+
+            if (nFaces > 0)
             {
-                Pout<< "  <<Writing " << faces.size()
+                Info<< "  <<Writing " << nFaces
                     << " skew faces to set " << faces.name() << endl;
                 faces.write();
             }
@@ -160,25 +180,29 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
         {
             //noFailedChecks++;
 
-            if (points.size() > 0)
+            label nPoints = returnReduce(points.size(), sumOp<label>());
+
+            if (nPoints > 0)
             {
-                Pout<< "  <<Writing " << points.size()
+                Info<< "  <<Writing " << nPoints
                     << " points on short edges to set " << points.name()
                     << endl;
                 points.write();
             }
         }
 
-        label nEdgeClose = points.size();
+        label nEdgeClose = returnReduce(points.size(), sumOp<label>());
 
         if (mesh.checkPointNearness(false, minDistSqr, &points))
         {
             //noFailedChecks++;
 
-            if (points.size() > nEdgeClose)
+            label nPoints = returnReduce(points.size(), sumOp<label>());
+
+            if (nPoints > nEdgeClose)
             {
                 pointSet nearPoints(mesh, "nearPoints", points);
-                Pout<< "  <<Writing " << nearPoints.size()
+                Info<< "  <<Writing " << nPoints
                     << " near (closer than " << Foam::sqrt(minDistSqr)
                     << " apart) points to set " << nearPoints.name() << endl;
                 nearPoints.write();
@@ -193,9 +217,11 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
         {
             //noFailedChecks++;
 
-            if (faces.size() > 0)
+            label nFaces = returnReduce(faces.size(), sumOp<label>());
+
+            if (nFaces > 0)
             {
-                Pout<< "  <<Writing " << faces.size()
+                Info<< "  <<Writing " << nFaces
                     << " faces with concave angles to set " << faces.name()
                     << endl;
                 faces.write();
@@ -210,9 +236,11 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
         {
             //noFailedChecks++;
 
-            if (faces.size() > 0)
+            label nFaces = returnReduce(faces.size(), sumOp<label>());
+
+            if (nFaces > 0)
             {
-                Pout<< "  <<Writing " << faces.size()
+                Info<< "  <<Writing " << nFaces
                     << " warped faces to set " << faces.name() << endl;
                 faces.write();
             }
@@ -226,7 +254,9 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
         {
             noFailedChecks++;
 
-            Pout<< "  <<Writing " << cells.size()
+            label nCells = returnReduce(cells.size(), sumOp<label>());
+
+            Info<< "  <<Writing " << nCells
                 << " under-determined cells to set " << cells.name() << endl;
             cells.write();
         }
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
index 85dc0bae62..ac924136e1 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
@@ -16,7 +16,7 @@ Foam::label Foam::checkTopology
 {
     label noFailedChecks = 0;
 
-    Pout<< "Checking topology..." << endl;
+    Info<< "Checking topology..." << endl;
 
     // Check if the boundary definition is unique
     mesh.boundaryMesh().checkDefinition(true);
@@ -30,7 +30,9 @@ Foam::label Foam::checkTopology
         {
             noFailedChecks++;
 
-            Pout<< "  <<Writing " << points.size()
+            label nPoints = returnReduce(points.size(), sumOp<label>());
+
+            Info<< "  <<Writing " << nPoints
                 << " unused points to set " << points.name() << endl;
             points.write();
         }
@@ -42,9 +44,12 @@ Foam::label Foam::checkTopology
         {
             noFailedChecks++;
         }
-        if (faces.size() > 0)
+
+        label nFaces = returnReduce(faces.size(), sumOp<label>());
+
+        if (nFaces > 0)
         {
-            Pout<< "  <<Writing " << faces.size()
+            Info<< "  <<Writing " << nFaces
                 << " unordered faces to set " << faces.name() << endl;
             faces.write();
         }
@@ -57,7 +62,9 @@ Foam::label Foam::checkTopology
         {
             noFailedChecks++;
 
-            Pout<< "  <<Writing " << cells.size()
+            label nCells = returnReduce(cells.size(), sumOp<label>());
+
+            Info<< "  <<Writing " << nCells
                 << " cells with over used edges to set " << cells.name()
                 << endl;
             cells.write();
@@ -70,7 +77,9 @@ Foam::label Foam::checkTopology
         {
             noFailedChecks++;
 
-            Pout<< "  <<Writing " << faces.size()
+            label nFaces = returnReduce(faces.size(), sumOp<label>());
+
+            Info<< "  <<Writing " << nFaces
                 << " faces with out-of-range or duplicate vertices to set "
                 << faces.name() << endl;
             faces.write();
@@ -84,7 +93,9 @@ Foam::label Foam::checkTopology
         {
             noFailedChecks++;
 
-            Pout<< "  <<Writing " << faces.size()
+            label nFaces = returnReduce(faces.size(), sumOp<label>());
+
+            Info<< "  <<Writing " << nFaces
                 << " faces with incorrect edges to set " << faces.name()
                 << endl;
             faces.write();
diff --git a/applications/utilities/mesh/manipulation/checkMesh/printMeshStats.C b/applications/utilities/mesh/manipulation/checkMesh/printMeshStats.C
index 9c608b689d..9a2c9396af 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/printMeshStats.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/printMeshStats.C
@@ -12,43 +12,72 @@
 
 void Foam::printMeshStats(const polyMesh& mesh, const bool allTopology)
 {
-    Pout<< "Mesh stats" << nl
-        << "    points:           " << mesh.points().size() << nl;
+    Info<< "Mesh stats" << nl
+        << "    points:           "
+        << returnReduce(mesh.points().size(), sumOp<label>()) << nl;
 
-    if (mesh.nInternalPoints() != -1)
+    label nInternalPoints = returnReduce
+    (
+        mesh.nInternalPoints(),
+        sumOp<label>()
+    );
+
+    if (nInternalPoints != -Pstream::nProcs())
     {
-        Pout<< "    internal points:  " << mesh.nInternalPoints() << nl;
+        Info<< "    internal points:  " << nInternalPoints << nl;
+
+        if (returnReduce(mesh.nInternalPoints(), minOp<label>()) == -1)
+        {
+            WarningIn("Foam::printMeshStats(const polyMesh&, const bool)")
+                << "Some processors have their points sorted into internal"
+                << " and external and some do not." << endl
+                << "This can cause problems later on." << endl;
+        }
     }
 
-    if (allTopology && mesh.nInternalPoints() != -1)
+    if (allTopology && nInternalPoints != -Pstream::nProcs())
     {
-        Pout<< "    edges:            " << mesh.nEdges() << nl
-            << "    internal edges:   " << mesh.nInternalEdges() << nl
+        label nEdges = returnReduce(mesh.nEdges(), sumOp<label>());
+        label nInternalEdges = returnReduce
+        (
+            mesh.nInternalEdges(),
+            sumOp<label>()
+        );
+        label nInternal1Edges = returnReduce
+        (
+            mesh.nInternal1Edges(),
+            sumOp<label>()
+        );
+        label nInternal0Edges = returnReduce
+        (
+            mesh.nInternal0Edges(),
+            sumOp<label>()
+        );
+
+        Info<< "    edges:            " << nEdges << nl
+            << "    internal edges:   " << nInternalEdges << nl
             << "    internal edges using one boundary point:   "
-            << mesh.nInternal1Edges()-mesh.nInternal0Edges() << nl
+            << nInternal1Edges-nInternal0Edges << nl
             << "    internal edges using two boundary points:  "
-            << mesh.nInternalEdges()-mesh.nInternal1Edges() << nl;
+            << nInternalEdges-nInternal1Edges << nl;
     }
 
-    Pout<< "    faces:            " << mesh.faces().size() << nl
-        << "    internal faces:   " << mesh.faceNeighbour().size() << nl
-        << "    cells:            " << mesh.cells().size() << nl
-        << "    boundary patches: " << mesh.boundaryMesh().size() << nl
-        << "    point zones:      " << mesh.pointZones().size() << nl
-        << "    face zones:       " << mesh.faceZones().size() << nl
-        << "    cell zones:       " << mesh.cellZones().size() << nl
+    Info<< "    faces:            "
+        << returnReduce(mesh.faces().size(), sumOp<label>()) << nl
+        << "    internal faces:   "
+        << returnReduce(mesh.faceNeighbour().size(), sumOp<label>()) << nl
+        << "    cells:            "
+        << returnReduce(mesh.cells().size(), sumOp<label>()) << nl
+        << "    boundary patches: "
+        << returnReduce(mesh.boundaryMesh().size(), sumOp<label>()) << nl
+        << "    point zones:      "
+        << returnReduce(mesh.pointZones().size(), sumOp<label>()) << nl
+        << "    face zones:       "
+        << returnReduce(mesh.faceZones().size(), sumOp<label>()) << nl
+        << "    cell zones:       "
+        << returnReduce(mesh.cellZones().size(), sumOp<label>()) << nl
         << endl;
 
-    if (Pstream::parRun())
-    {
-        const globalMeshData& parData = mesh.globalData();
-
-        Info<< "Overall stats" << nl
-            << "    points:   " << parData.nTotalPoints() << nl
-            << "    faces:    " << parData.nTotalFaces() << nl
-            << "    cells:    " << parData.nTotalCells() << nl
-            << endl;
-    }
 
     // Construct shape recognizers
     hexMatcher hex;
@@ -99,7 +128,15 @@ void Foam::printMeshStats(const polyMesh& mesh, const bool allTopology)
         }
     }
 
-    Pout<< "Number of cells of each type: " << nl
+    reduce(nHex,sumOp<label>());
+    reduce(nPrism,sumOp<label>()); 
+    reduce(nWedge,sumOp<label>());
+    reduce(nPyr,sumOp<label>());
+    reduce(nTetWedge,sumOp<label>());
+    reduce(nTet,sumOp<label>());
+    reduce(nUnknown,sumOp<label>());
+
+    Info<< "Overall number of cells of each type:" << nl
         << "    hexahedra:     " << nHex << nl
         << "    prisms:        " << nPrism << nl
         << "    wedges:        " << nWedge << nl
diff --git a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
index a74458fe9f..e7d7567a13 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
@@ -552,6 +552,8 @@ bool Foam::polyBoundaryMesh::checkDefinition(const bool report) const
         nextPatchStart += bm[patchI].size();
     }
 
+    reduce(boundaryError, orOp<bool>());
+
     if (boundaryError)
     {
         if (debug || report)
@@ -565,7 +567,7 @@ bool Foam::polyBoundaryMesh::checkDefinition(const bool report) const
     {
         if (debug || report)
         {
-            Pout << "    Boundary definition OK." << endl;
+            Info << "    Boundary definition OK." << endl;
         }
 
         return false;

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 06/10] 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;
 
 // ************************************************************************* //

From 3feae1f3f6955ae8cb9228c325e4db00de99acda Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Wed, 8 Oct 2008 11:41:38 +0100
Subject: [PATCH 07/10] extraneous copy

---
 .../manipulation/checkMesh.save/Make/files    |   6 -
 .../manipulation/checkMesh.save/Make/options  |   5 -
 .../checkMesh.save/checkGeometry.C            | 220 ----------------
 .../checkMesh.save/checkGeometry.H            |  13 -
 .../manipulation/checkMesh.save/checkMesh.C   | 152 -----------
 .../checkMesh.save/checkTopology.C            | 236 ------------------
 .../checkMesh.save/checkTopology.H            |   8 -
 .../checkMesh.save/printMeshStats.C           |  95 -------
 .../checkMesh.save/printMeshStats.H           |   6 -
 9 files changed, 741 deletions(-)
 delete mode 100644 applications/utilities/mesh/manipulation/checkMesh.save/Make/files
 delete mode 100644 applications/utilities/mesh/manipulation/checkMesh.save/Make/options
 delete mode 100644 applications/utilities/mesh/manipulation/checkMesh.save/checkGeometry.C
 delete mode 100644 applications/utilities/mesh/manipulation/checkMesh.save/checkGeometry.H
 delete mode 100644 applications/utilities/mesh/manipulation/checkMesh.save/checkMesh.C
 delete mode 100644 applications/utilities/mesh/manipulation/checkMesh.save/checkTopology.C
 delete mode 100644 applications/utilities/mesh/manipulation/checkMesh.save/checkTopology.H
 delete mode 100644 applications/utilities/mesh/manipulation/checkMesh.save/printMeshStats.C
 delete mode 100644 applications/utilities/mesh/manipulation/checkMesh.save/printMeshStats.H

diff --git a/applications/utilities/mesh/manipulation/checkMesh.save/Make/files b/applications/utilities/mesh/manipulation/checkMesh.save/Make/files
deleted file mode 100644
index f0b7c16694..0000000000
--- a/applications/utilities/mesh/manipulation/checkMesh.save/Make/files
+++ /dev/null
@@ -1,6 +0,0 @@
-printMeshStats.C
-checkTopology.C
-checkGeometry.C
-checkMesh.C
-
-EXE = $(FOAM_APPBIN)/checkMesh
diff --git a/applications/utilities/mesh/manipulation/checkMesh.save/Make/options b/applications/utilities/mesh/manipulation/checkMesh.save/Make/options
deleted file mode 100644
index 54c035b8f5..0000000000
--- a/applications/utilities/mesh/manipulation/checkMesh.save/Make/options
+++ /dev/null
@@ -1,5 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/meshTools/lnInclude
-
-EXE_LIBS = \
-    -lmeshTools
diff --git a/applications/utilities/mesh/manipulation/checkMesh.save/checkGeometry.C b/applications/utilities/mesh/manipulation/checkMesh.save/checkGeometry.C
deleted file mode 100644
index 9eb4ef63ad..0000000000
--- a/applications/utilities/mesh/manipulation/checkMesh.save/checkGeometry.C
+++ /dev/null
@@ -1,220 +0,0 @@
-#include "checkGeometry.H"
-#include "polyMesh.H"
-#include "globalMeshData.H"
-#include "cellSet.H"
-#include "faceSet.H"
-#include "pointSet.H"
-
-Foam::label Foam::checkGeometry
-(
-    const polyMesh& mesh,
-    bool checkPointNearness,
-    bool checkCellDeterminant
-)
-{
-    label noFailedChecks = 0;
-
-    Info<< "\nChecking geometry..." << endl;
-
-    boundBox bb(mesh.points());
-
-    Pout<< "    Domain bounding box: "
-           << bb.min() << " " << bb.max() << endl;
-
-    // Get a small relative length from the bounding box
-    const boundBox& globalBb = mesh.globalData().bb();
-
-    if (Pstream::parRun())
-    {
-        Info<< "    Overall domain bounding box: "
-            << globalBb.min() << " " << globalBb.max() << endl;
-    }
-
-
-    // Min length
-    scalar minDistSqr = magSqr(1e-6*(globalBb.max() - globalBb.min()));
-
-    
-    if (mesh.checkClosedBoundary(true)) noFailedChecks++;
-
-    {
-        cellSet cells(mesh, "nonClosedCells", mesh.nCells()/100+1);
-        cellSet aspectCells(mesh, "highAspectRatioCells", mesh.nCells()/100+1);
-        if (mesh.checkClosedCells(true, &cells, &aspectCells))
-        {
-            noFailedChecks++;
-
-            if (cells.size() > 0)
-            {
-                Pout<< "  <<Writing " << cells.size()
-                    << " non closed cells to set " << cells.name() << endl;
-                cells.write();
-            }
-        }
-        if (aspectCells.size() > 0)
-        {
-            Pout<< "  <<Writing " << aspectCells.size()
-                << " cells with high aspect ratio to set "
-                << aspectCells.name() << endl;
-            aspectCells.write();
-        }
-    }
-
-    {
-        faceSet faces(mesh, "zeroAreaFaces", mesh.nFaces()/100 + 1);
-        if (mesh.checkFaceAreas(true, &faces))
-        {
-            noFailedChecks++;
-
-            if (faces.size() > 0)
-            {
-                Pout<< "  <<Writing " << faces.size()
-                    << " zero area faces to set " << faces.name() << endl;
-                faces.write();
-            }
-        }
-    }
-
-    {
-        cellSet cells(mesh, "zeroVolumeCells", mesh.nCells()/100 + 1);
-        if (mesh.checkCellVolumes(true, &cells))
-        {
-            noFailedChecks++;
-
-            if (cells.size() > 0)
-            {
-                Pout<< "  <<Writing " << cells.size()
-                    << " zero volume cells to set " << cells.name() << endl;
-                cells.write();
-            }
-        }
-    }
-
-    {
-        faceSet faces(mesh, "nonOrthoFaces", mesh.nFaces()/100 + 1);
-        if (mesh.checkFaceOrthogonality(true, &faces))
-        {
-            noFailedChecks++;
-        }
-
-        if (faces.size() > 0)
-        {
-            Pout<< "  <<Writing " << faces.size()
-                << " non-orthogonal faces to set " << faces.name() << endl;
-            faces.write();
-        }
-    }
-
-
-    {
-        faceSet faces(mesh, "wrongOrientedFaces", mesh.nFaces()/100 + 1);
-        if (mesh.checkFacePyramids(true, -SMALL, &faces))
-        {
-            noFailedChecks++;
-
-            if (faces.size() > 0)
-            {
-                Pout<< "  <<Writing " << faces.size()
-                    << " faces with incorrect orientation to set "
-                    << faces.name() << endl;
-                faces.write();
-            }
-        }
-    }
-
-    {
-        faceSet faces(mesh, "skewFaces", mesh.nFaces()/100 + 1);
-        if (mesh.checkFaceSkewness(true, &faces))
-        {
-            noFailedChecks++;
-
-            if (faces.size() > 0)
-            {
-                Pout<< "  <<Writing " << faces.size()
-                    << " skew faces to set " << faces.name() << endl;
-                faces.write();
-            }
-        }
-    }
-
-    if (checkPointNearness)
-    {
-        // Note use of nPoints since don't want edge construction.
-        pointSet points(mesh, "shortEdges", mesh.nPoints()/1000 + 1);
-        if (mesh.checkEdgeLength(true, minDistSqr, &points))
-        {
-            //noFailedChecks++;
-
-            if (points.size() > 0)
-            {
-                Pout<< "  <<Writing " << points.size()
-                    << " points on short edges to set " << points.name()
-                    << endl;
-                points.write();
-            }
-        }
-
-        label nEdgeClose = points.size();
-
-        if (mesh.checkPointNearness(false, minDistSqr, &points))
-        {
-            //noFailedChecks++;
-
-            if (points.size() > nEdgeClose)
-            {
-                pointSet nearPoints(mesh, "nearPoints", points);
-                Pout<< "  <<Writing " << nearPoints.size()
-                    << " near (closer than " << Foam::sqrt(minDistSqr)
-                    << " apart) points to set " << nearPoints.name() << endl;
-                nearPoints.write();
-            }
-        }
-    }
-
-    {
-        faceSet faces(mesh, "concaveFaces", mesh.nFaces()/100 + 1);
-        if (mesh.checkFaceAngles(true, 10, &faces))
-        {
-            //noFailedChecks++;
-
-            if (faces.size() > 0)
-            {
-                Pout<< "  <<Writing " << faces.size()
-                    << " faces with concave angles to set " << faces.name()
-                    << endl;
-                faces.write();
-            }
-        }
-    }
-
-    {
-        faceSet faces(mesh, "warpedFaces", mesh.nFaces()/100 + 1);
-        if (mesh.checkFaceFlatness(true, 0.8, &faces))
-        {
-            //noFailedChecks++;
-
-            if (faces.size() > 0)
-            {
-                Pout<< "  <<Writing " << faces.size()
-                    << " warped faces to set " << faces.name() << endl;
-                faces.write();
-            }
-        }
-    }
-
-    if (checkCellDeterminant)
-    {
-        cellSet cells(mesh, "underdeterminedCells", mesh.nCells()/100);
-        if (mesh.checkCellDeterminant(true, &cells))
-        {
-            noFailedChecks++;
-
-            Pout<< "  <<Writing " << cells.size()
-                << " under-determines cells to set " << cells.name() << endl;
-            cells.write();
-        }
-    }
-    
-
-    return noFailedChecks;
-}
diff --git a/applications/utilities/mesh/manipulation/checkMesh.save/checkGeometry.H b/applications/utilities/mesh/manipulation/checkMesh.save/checkGeometry.H
deleted file mode 100644
index 791b007820..0000000000
--- a/applications/utilities/mesh/manipulation/checkMesh.save/checkGeometry.H
+++ /dev/null
@@ -1,13 +0,0 @@
-#include "label.H"
-
-namespace Foam
-{
-    class polyMesh;
-
-    label checkGeometry
-    (
-        const polyMesh& mesh,
-        bool checkPointNearness,
-        bool checkCellDeterminant
-    );
-}
diff --git a/applications/utilities/mesh/manipulation/checkMesh.save/checkMesh.C b/applications/utilities/mesh/manipulation/checkMesh.save/checkMesh.C
deleted file mode 100644
index ec598b4c21..0000000000
--- a/applications/utilities/mesh/manipulation/checkMesh.save/checkMesh.C
+++ /dev/null
@@ -1,152 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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
-
-Application
-    checkMesh
-
-Description
-    Checks validity of a mesh
-
-\*---------------------------------------------------------------------------*/
-
-#include "argList.H"
-#include "Time.H"
-#include "polyMesh.H"
-#include "globalMeshData.H"
-
-#include "printMeshStats.H"
-#include "checkTopology.H"
-#include "checkGeometry.H"
-
-using namespace Foam;
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-
-#   include "addTimeOptionsNoConstant.H"
-
-    argList::validOptions.insert("fullTopology", "");
-    argList::validOptions.insert("pointNearness", "");
-    argList::validOptions.insert("cellDeterminant", "");
-
-#   include "setRootCase.H"
-#   include "createTime.H"
-
-    // Get times list
-    instantList Times = runTime.times();
-
-#   include "checkTimeOptionsNoConstant.H"
-
-    runTime.setTime(Times[startTime], startTime);
-
-#   include "createPolyMesh.H"
-
-    bool firstCheck = true;
-
-    for (label i=startTime; i<endTime; i++)
-    {
-        runTime.setTime(Times[i], i);
-
-        polyMesh::readUpdateState state = mesh.readUpdate();
-
-        if
-        (
-            firstCheck
-         || state == polyMesh::TOPO_CHANGE
-         || state == polyMesh::TOPO_PATCH_CHANGE
-        )
-        {
-            firstCheck = false;
-
-            Info<< "Time = " << runTime.timeName() << nl << endl;
-
-            // Clear mesh before checking
-            mesh.clearOut();
-
-            // Reconstruct globalMeshData
-            mesh.globalData();
-
-            printMeshStats(mesh);
-
-            label noFailedChecks = 0;
-
-            noFailedChecks += checkTopology
-            (
-                mesh,
-                args.options().found("fullTopology")
-            );
-
-            noFailedChecks += checkGeometry
-            (
-                mesh,
-                args.options().found("pointNearness"),
-                args.options().found("cellDeterminant")
-            );
-
-            reduce(noFailedChecks, sumOp<label>());
-
-            if (noFailedChecks == 0)
-            {
-                Info<< "\nMesh OK."
-                    << nl << endl;
-            }
-            else
-            {
-                Info<< "\nFailed " << noFailedChecks << " mesh checks."
-                    << nl << endl;
-            }
-        }
-        else if (state == polyMesh::POINTS_MOVED)
-        {
-            label noFailedChecks = checkGeometry
-            (
-                mesh,
-                args.options().found("pointNearness"),
-                args.options().found("cellDeterminant")
-            );
-
-            reduce(noFailedChecks, sumOp<label>());
-
-            if (noFailedChecks == 0)
-            {
-                Info << "\nMesh OK."
-                    << nl << endl;
-            }
-            else
-            {
-                Info<< "\nFailed " << noFailedChecks << " mesh checks."
-                    << nl << endl;
-            }
-        }
-    }
-
-    Info<< "End\n" << endl;
-
-    return(0);
-}
-
-
-// ************************************************************************* //
diff --git a/applications/utilities/mesh/manipulation/checkMesh.save/checkTopology.C b/applications/utilities/mesh/manipulation/checkMesh.save/checkTopology.C
deleted file mode 100644
index 885c5f36aa..0000000000
--- a/applications/utilities/mesh/manipulation/checkMesh.save/checkTopology.C
+++ /dev/null
@@ -1,236 +0,0 @@
-#include "checkTopology.H"
-#include "polyMesh.H"
-#include "Time.H"
-#include "regionSplit.H"
-#include "cellSet.H"
-#include "faceSet.H"
-#include "pointSet.H"
-#include "IOmanip.H"
-
-Foam::label Foam::checkTopology(const polyMesh& mesh, bool fullTopology)
-{
-    label noFailedChecks = 0;
-
-    Pout<< "Checking topology..." << endl;
-
-    // Check if the boundary definition is unique
-    mesh.boundaryMesh().checkDefinition(true);
-
-    // Check if the boundary processor patches are correct
-    mesh.boundaryMesh().checkParallelSync(true);
-
-    {
-        pointSet points(mesh, "unusedPoints", mesh.nPoints()/100);
-        if (mesh.checkPoints(true, &points))
-        {
-            noFailedChecks++;
-
-            Pout<< "  <<Writing " << points.size()
-                << " unused points to set " << points.name() << endl;
-            points.write();
-        }
-    }
-
-    {
-        faceSet faces(mesh, "upperTriangularFace", mesh.nFaces()/100);
-        if (mesh.checkUpperTriangular(true, &faces))
-        {
-            noFailedChecks++;
-
-            Pout<< "  <<Writing " << faces.size()
-                << " unordered faces to set " << faces.name() << endl;
-            faces.write();
-        }
-    }
-
-    {
-        cellSet cells(mesh, "zipUpCells", mesh.nCells()/100);
-        if (mesh.checkCellsZipUp(true, &cells))
-        {
-            noFailedChecks++;
-
-            Pout<< "  <<Writing " << cells.size()
-                << " cells with over used edges to set " << cells.name()
-                << endl;
-            cells.write();
-        }
-    }
-
-    {
-        faceSet faces(mesh, "outOfRangeFaces", mesh.nFaces()/100);
-        if (mesh.checkFaceVertices(true, &faces))
-        {
-            noFailedChecks++;
-
-            Pout<< "  <<Writing " << faces.size()
-                << " faces with out-of-range vertices to set " << faces.name()
-                << endl;
-            faces.write();
-        }
-    }
-
-    {
-        faceSet faces(mesh, "edgeFaces", mesh.nFaces()/100);
-        if (mesh.checkFaceFaces(true, &faces))
-        {
-            noFailedChecks++;
-
-            Pout<< "  <<Writing " << faces.size()
-                << " faces with incorrect edges to set " << faces.name()
-                << endl;
-            faces.write();
-        }
-    }
-
-    {
-        regionSplit rs(mesh);
-
-        if (rs.nRegions() == 1)
-        {
-            Info<< "    Number of regions: " << rs.nRegions() << " (OK)."
-                << endl;
-        
-        }
-        else
-        {
-            Info<< "   *Number of regions: " << rs.nRegions() << endl;
-
-            Info<< "    The mesh has multiple regions which are not connected "
-                   "by any face." << endl
-                << "  <<Writing region information to "
-                << mesh.time().timeName()/"cellToRegion"
-                << endl;
-
-            labelIOList ctr
-            (
-                IOobject
-                (
-                    "cellToRegion",
-                    mesh.time().timeName(),
-                    mesh,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
-                rs
-            );
-            ctr.write();
-        }
-    }
-
-    if (!Pstream::parRun())
-    {
-        Pout<< "\nChecking patch topology for multiply connected surfaces ..."
-            << endl;
-
-        const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-        // Non-manifold points
-        pointSet points
-        (
-            mesh,
-            "nonManifoldPoints",
-            mesh.nPoints()/100
-        );
-
-        Pout.setf(ios_base::left);
-
-        Pout<< "    "
-            << setw(20) << "Patch"
-            << setw(9) << "Faces"
-            << setw(9) << "Points"
-            << " Surface" << endl;
-
-        forAll(patches, patchI)
-        {
-            const polyPatch& pp = patches[patchI];
-
-            primitivePatch::surfaceTopo pTyp = pp.surfaceType();
-
-            if (pp.size() == 0)
-            {
-                Pout<< "    "
-                    << setw(20) << pp.name()
-                    << setw(9) << pp.size()
-                    << setw(9) << pp.nPoints()
-                    << " ok (empty)" << endl;
-            }
-            else if (pTyp == primitivePatch::MANIFOLD)
-            {
-                if (pp.checkPointManifold(true, &points))
-                {
-                    Pout<< "    " 
-                        << setw(20) << pp.name()
-                        << setw(9) << pp.size()
-                        << setw(9) << pp.nPoints()
-                        << " multiply connected (shared point)" << endl;
-                }
-                else
-                {
-                    Pout<< "    "
-                        << setw(20) << pp.name()
-                        << setw(9) << pp.size()
-                        << setw(9) << pp.nPoints()
-                        << " ok (closed singly connected surface)" << endl;
-                }
-
-                // Add points on non-manifold edges to make set complete
-                pp.checkTopology(false, &points);
-            }
-            else
-            {
-                pp.checkTopology(false, &points);
-
-                if (pTyp == primitivePatch::OPEN)
-                {
-                    Pout<< "    "
-                        << setw(20) << pp.name()
-                        << setw(9) << pp.size()
-                        << setw(9) << pp.nPoints()
-                        << " ok (not multiply connected)" << endl;
-                }
-                else
-                {
-                    Pout<< "    "
-                        << setw(20) << pp.name()
-                        << setw(9) << pp.size()
-                        << setw(9) << pp.nPoints()
-                        << " multiply connected surface (shared edge)"
-                        << endl;
-                }
-            }
-        }
-
-        if (points.size() > 0)
-        {
-            Pout<< "  <<Writing " << points.size()
-                << " conflicting points to set "
-                << points.name() << endl;
-
-            points.write();
-        }
-
-        //Pout.setf(ios_base::right);
-    }
-
-    // Force creation of all addressing if requested.
-    // Errors will be reported as required
-    if (fullTopology)
-    {
-        mesh.cells();
-        mesh.faces();
-        mesh.edges();
-        mesh.points();
-        mesh.faceOwner();
-        mesh.faceNeighbour();
-        mesh.cellCells();
-        mesh.edgeCells();
-        mesh.pointCells();
-        mesh.edgeFaces();
-        mesh.pointFaces();
-        mesh.cellEdges();
-        mesh.faceEdges();
-        mesh.pointEdges();
-    }
-
-    return noFailedChecks;
-}
diff --git a/applications/utilities/mesh/manipulation/checkMesh.save/checkTopology.H b/applications/utilities/mesh/manipulation/checkMesh.save/checkTopology.H
deleted file mode 100644
index 6db7e0d3c5..0000000000
--- a/applications/utilities/mesh/manipulation/checkMesh.save/checkTopology.H
+++ /dev/null
@@ -1,8 +0,0 @@
-#include "label.H"
-
-namespace Foam
-{
-    class polyMesh;
-
-    label checkTopology(const polyMesh& mesh, bool fullTopology);
-}
diff --git a/applications/utilities/mesh/manipulation/checkMesh.save/printMeshStats.C b/applications/utilities/mesh/manipulation/checkMesh.save/printMeshStats.C
deleted file mode 100644
index 30da859db6..0000000000
--- a/applications/utilities/mesh/manipulation/checkMesh.save/printMeshStats.C
+++ /dev/null
@@ -1,95 +0,0 @@
-#include "printMeshStats.H"
-#include "polyMesh.H"
-#include "globalMeshData.H"
-
-#include "hexMatcher.H"
-#include "wedgeMatcher.H"
-#include "prismMatcher.H"
-#include "pyrMatcher.H"
-#include "tetWedgeMatcher.H"
-#include "tetMatcher.H"
-
-
-void Foam::printMeshStats(const polyMesh& mesh)
-{
-    Pout<< "Mesh stats" << nl
-        << "    points:           " << mesh.points().size() << nl
-        << "    faces:            " << mesh.faces().size() << nl
-        << "    internal faces:   " << mesh.faceNeighbour().size() << nl
-        << "    cells:            " << mesh.cells().size() << nl
-        << "    boundary patches: " << mesh.boundaryMesh().size() << nl
-        << "    point zones:      " << mesh.pointZones().size() << nl
-        << "    face zones:       " << mesh.faceZones().size() << nl
-        << "    cell zones:       " << mesh.cellZones().size() << nl
-        << endl;
-
-    if (Pstream::parRun())
-    {
-        const globalMeshData& parData = mesh.globalData();
-
-        Info<< "Overall stats" << nl
-            << "    points:   " << parData.nTotalPoints() << nl
-            << "    faces:    " << parData.nTotalFaces() << nl
-            << "    cells:    " << parData.nTotalCells() << nl
-            << endl;
-    }
-
-    // Construct shape recognizers
-    hexMatcher hex;
-    prismMatcher prism;
-    wedgeMatcher wedge;
-    pyrMatcher pyr;
-    tetWedgeMatcher tetWedge;
-    tetMatcher tet;
-
-    // Counters for different cell types
-    label nHex = 0;
-    label nWedge = 0;
-    label nPrism = 0;
-    label nPyr = 0;
-    label nTet = 0;
-    label nTetWedge = 0;
-    label nUnknown = 0;
-
-    for(label cellI = 0; cellI < mesh.nCells(); cellI++)
-    {
-        if (hex.isA(mesh, cellI))
-        {
-            nHex++;
-        }
-        else if (tet.isA(mesh, cellI))
-        {
-            nTet++;
-        }
-        else if (pyr.isA(mesh, cellI))
-        {
-            nPyr++;
-        }
-        else if (prism.isA(mesh, cellI))
-        {
-            nPrism++;
-        }
-        else if (wedge.isA(mesh, cellI))
-        {
-            nWedge++;
-        }
-        else if (tetWedge.isA(mesh, cellI))
-        {
-            nTetWedge++;
-        }
-        else
-        {   
-            nUnknown++;
-        }
-    }
-
-    Pout<< "Number of cells of each type: " << nl
-        << "    hexahedra:     " << nHex << nl
-        << "    prisms:        " << nPrism << nl
-        << "    wedges:        " << nWedge << nl
-        << "    pyramids:      " << nPyr << nl
-        << "    tet wedges:    " << nTetWedge << nl
-        << "    tetrahedra:    " << nTet << nl
-        << "    polyhedra:     " << nUnknown
-        << nl << endl;
-}
diff --git a/applications/utilities/mesh/manipulation/checkMesh.save/printMeshStats.H b/applications/utilities/mesh/manipulation/checkMesh.save/printMeshStats.H
deleted file mode 100644
index 2da1201188..0000000000
--- a/applications/utilities/mesh/manipulation/checkMesh.save/printMeshStats.H
+++ /dev/null
@@ -1,6 +0,0 @@
-namespace Foam
-{
-    class polyMesh;
-
-    void printMeshStats(const polyMesh& mesh);
-}

From 12bec33d048158c62365dfe570002e7728fe4f68 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Wed, 8 Oct 2008 11:42:23 +0100
Subject: [PATCH 08/10] translation vector on cyclicPolyPatch; adapted
 createPatch

---
 .../manipulation/createPatch/createPatch.C    | 375 ++++++++++++++++--
 .../manipulation/createPatch/createPatchDict  |  63 ++-
 .../basic/coupled/coupledPolyPatch.C          |  10 +-
 .../constraint/cyclic/cyclicPolyPatch.C       | 173 +++++---
 .../constraint/cyclic/cyclicPolyPatch.H       | 139 ++++---
 5 files changed, 601 insertions(+), 159 deletions(-)

diff --git a/applications/utilities/mesh/manipulation/createPatch/createPatch.C b/applications/utilities/mesh/manipulation/createPatch/createPatch.C
index a449661982..01df2c2f28 100644
--- a/applications/utilities/mesh/manipulation/createPatch/createPatch.C
+++ b/applications/utilities/mesh/manipulation/createPatch/createPatch.C
@@ -34,6 +34,7 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
+#include "cyclicPolyPatch.H"
 #include "syncTools.H"
 #include "argList.H"
 #include "polyMesh.H"
@@ -256,27 +257,6 @@ void dumpCyclicMatch(const fileName& prefix, const polyMesh& mesh)
                 );
             }
 
-//            cycPatch.writeOBJ
-//            (
-//                prefix+cycPatch.name()+"_half0.obj",
-//                SubList<face>
-//                (
-//                    cycPatch,
-//                    halfSize
-//                ),
-//                cycPatch.points()
-//            );
-//            cycPatch.writeOBJ
-//            (
-//                prefix+cycPatch.name()+"_half1.obj",
-//                SubList<face>
-//                (
-//                    cycPatch,
-//                    halfSize,
-//                    halfSize
-//                ),
-//                cycPatch.points()
-//            );
 
             // Lines between corresponding face centres
             OFstream str(prefix+cycPatch.name()+"_match.obj");
@@ -289,7 +269,8 @@ void dumpCyclicMatch(const fileName& prefix, const polyMesh& mesh)
                 vertI++;
 
                 label nbrFaceI = halfSize + faceI;
-                const point& fc1 = mesh.faceCentres()[cycPatch.start()+nbrFaceI];
+                const point& fc1 =
+                    mesh.faceCentres()[cycPatch.start()+nbrFaceI];
                 meshTools::writeOBJ(str, fc1);
                 vertI++;
 
@@ -300,6 +281,247 @@ void dumpCyclicMatch(const fileName& prefix, const polyMesh& mesh)
 }
 
 
+void separateList
+(
+    const vectorField& separation,
+    UList<vector>& field
+)
+{
+    if (separation.size() == 1)
+    {
+        // Single value for all.
+
+        forAll(field, i)
+        {
+            field[i] += separation[0];
+        }
+    }
+    else if (separation.size() == field.size())
+    {
+        forAll(field, i)
+        {
+            field[i] += separation[i];
+        }
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "separateList(const vectorField&, UList<vector>&)"
+        )   << "Sizes of field and transformation not equal. field:"
+            << field.size() << " transformation:" << separation.size()
+            << abort(FatalError);
+    }
+}
+
+
+// Synchronise points on both sides of coupled boundaries.
+template <class CombineOp>
+void syncPoints
+(
+    const polyMesh& mesh,
+    pointField& points,
+    const CombineOp& cop,
+    const point& nullValue
+)
+{
+    if (points.size() != mesh.nPoints())
+    {
+        FatalErrorIn
+        (
+            "syncPoints"
+            "(const polyMesh&, pointField&, const CombineOp&, const point&)"
+        )   << "Number of values " << points.size()
+            << " is not equal to the number of points in the mesh "
+            << mesh.nPoints() << abort(FatalError);
+    }
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    // Is there any coupled patch with transformation?
+    bool hasTransformation = false;
+
+    if (Pstream::parRun())
+    {
+        // Send
+
+        forAll(patches, patchI)
+        {
+            const polyPatch& pp = patches[patchI];
+
+            if
+            (
+                isA<processorPolyPatch>(pp)
+             && pp.nPoints() > 0
+             && refCast<const processorPolyPatch>(pp).owner()
+            )
+            {
+                const processorPolyPatch& procPatch =
+                    refCast<const processorPolyPatch>(pp);
+
+                // Get data per patchPoint in neighbouring point numbers.
+                pointField patchInfo(procPatch.nPoints(), nullValue);
+
+                const labelList& meshPts = procPatch.meshPoints();
+                const labelList& nbrPts = procPatch.neighbPoints();
+
+                forAll(nbrPts, pointI)
+                {
+                    label nbrPointI = nbrPts[pointI];
+                    if (nbrPointI >= 0 && nbrPointI < patchInfo.size())
+                    {
+                        patchInfo[nbrPointI] = points[meshPts[pointI]];
+                    }
+                }
+
+                OPstream toNbr(Pstream::blocking, procPatch.neighbProcNo());
+                toNbr << patchInfo;
+            }
+        }
+
+
+        // Receive and set.
+
+        forAll(patches, patchI)
+        {
+            const polyPatch& pp = patches[patchI];
+
+            if
+            (
+                isA<processorPolyPatch>(pp)
+             && pp.nPoints() > 0
+             && !refCast<const processorPolyPatch>(pp).owner()
+            )
+            {
+                const processorPolyPatch& procPatch =
+                    refCast<const processorPolyPatch>(pp);
+
+                pointField nbrPatchInfo(procPatch.nPoints());
+                {
+                    // We do not know the number of points on the other side
+                    // so cannot use Pstream::read.
+                    IPstream fromNbr
+                    (
+                        Pstream::blocking,
+                        procPatch.neighbProcNo()
+                    );
+                    fromNbr >> nbrPatchInfo;
+                }
+                // Null any value which is not on neighbouring processor
+                nbrPatchInfo.setSize(procPatch.nPoints(), nullValue);
+
+                if (!procPatch.parallel())
+                {
+                    hasTransformation = true;
+                    transformList(procPatch.forwardT(), nbrPatchInfo);
+                }
+                else if (procPatch.separated())
+                {
+                    hasTransformation = true;
+                    separateList(-procPatch.separation(), nbrPatchInfo);
+                }
+
+                const labelList& meshPts = procPatch.meshPoints();
+
+                forAll(meshPts, pointI)
+                {
+                    label meshPointI = meshPts[pointI];
+                    points[meshPointI] = nbrPatchInfo[pointI];
+                }
+            }
+        }
+    }
+
+    // Do the cyclics.
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (isA<cyclicPolyPatch>(pp))
+        {
+            const cyclicPolyPatch& cycPatch =
+                refCast<const cyclicPolyPatch>(pp);
+
+            const edgeList& coupledPoints = cycPatch.coupledPoints();
+            const labelList& meshPts = cycPatch.meshPoints();
+
+            pointField half0Values(coupledPoints.size());
+
+            forAll(coupledPoints, i)
+            {
+                const edge& e = coupledPoints[i];
+                label point0 = meshPts[e[0]];
+                half0Values[i] = points[point0];
+            }
+
+            if (!cycPatch.parallel())
+            {
+                hasTransformation = true;
+                transformList(cycPatch.reverseT(), half0Values);
+            }
+            else if (cycPatch.separated())
+            {
+                hasTransformation = true;
+                const vectorField& v = cycPatch.coupledPolyPatch::separation();
+                separateList(v, half0Values);
+            }
+
+            forAll(coupledPoints, i)
+            {
+                const edge& e = coupledPoints[i];
+                label point1 = meshPts[e[1]];
+                points[point1] = half0Values[i];
+            }
+        }
+    }
+
+    //- Note: hasTransformation is only used for warning messages so
+    //  reduction not strictly nessecary.
+    //reduce(hasTransformation, orOp<bool>());
+
+    // Synchronize multiple shared points.
+    const globalMeshData& pd = mesh.globalData();
+
+    if (pd.nGlobalPoints() > 0)
+    {
+        if (hasTransformation)
+        {
+            WarningIn
+            (
+                "syncPoints"
+                "(const polyMesh&, pointField&, const CombineOp&, const point&)"
+            )   << "There are decomposed cyclics in this mesh with"
+                << " transformations." << endl
+                << "This is not supported. The result will be incorrect"
+                << endl;
+        }
+
+
+        // Values on shared points.
+        pointField sharedPts(pd.nGlobalPoints(), nullValue);
+
+        forAll(pd.sharedPointLabels(), i)
+        {
+            label meshPointI = pd.sharedPointLabels()[i];
+            // Fill my entries in the shared points
+            sharedPts[pd.sharedPointAddr()[i]] = points[meshPointI];
+        }
+
+        // Combine on master.
+        Pstream::listCombineGather(sharedPts, cop);
+        Pstream::listCombineScatter(sharedPts);
+
+        // Now we will all have the same information. Merge it back with
+        // my local information.
+        forAll(pd.sharedPointLabels(), i)
+        {
+            label meshPointI = pd.sharedPointLabels()[i];
+            points[meshPointI] = sharedPts[pd.sharedPointAddr()[i]];
+        }
+    }
+}
+
+
 // Main program:
 
 int main(int argc, char *argv[])
@@ -393,25 +615,28 @@ int main(int argc, char *argv[])
 
             label destPatchI = patches.findPatchID(patchName);
 
-            word patchType(dict.lookup("type"));
-
             if (destPatchI == -1)
             {
+                dictionary patchDict(dict.subDict("dictionary"));
+
                 destPatchI = allPatches.size();
 
                 Info<< "Adding new patch " << patchName
-                    << " of type " << patchType
-                    << " as patch " << destPatchI << endl;
+                    << " as patch " << destPatchI
+                    << " from " << patchDict << endl;
+
+                patchDict.remove("nFaces");
+                patchDict.add("nFaces", 0);
+                patchDict.remove("startFace");
+                patchDict.add("startFace", startFaceI);
 
                 // Add an empty patch.
                 allPatches.append
                 (
                     polyPatch::New
                     (
-                        patchType,
                         patchName,
-                        0,              // size
-                        startFaceI,     // start
+                        patchDict,
                         destPatchI,
                         patches
                     ).ptr()
@@ -557,16 +782,100 @@ int main(int argc, char *argv[])
     }
     else
     {
+        Info<< "Synchronising points." << nl << endl;
+
+        // This is a bit tricky. Both normal and position might be out and
+        // current separation also includes the normal
+        // ( separation_ = (nf&(Cr - Cf))*nf ).
+
+        // For processor patches:
+        // - disallow multiple separation/transformation. This basically
+        //   excludes decomposed cyclics. Use the (probably 0) separation
+        //   to align the points.
+        // For cyclic patches:
+        // - for separated ones use our own recalculated offset vector
+        // - for rotational ones use current one.
+
+        forAll(mesh.boundaryMesh(), patchI)
+        {
+            const polyPatch& pp = mesh.boundaryMesh()[patchI];
+
+            if (pp.size() && isA<coupledPolyPatch>(pp))
+            {
+                const coupledPolyPatch& cpp =
+                    refCast<const coupledPolyPatch>(pp);
+
+                if (cpp.separated())
+                {
+                    Info<< "On coupled patch " << pp.name()
+                        << " separation[0] was "
+                        << cpp.separation()[0] << endl;
+
+                    if (isA<cyclicPolyPatch>(pp))
+                    {
+                        const cyclicPolyPatch& cycpp =
+                            refCast<const cyclicPolyPatch>(pp);
+
+                        if (cycpp.transform() == cyclicPolyPatch::TRANSLATIONAL)
+                        {
+                            Info<< "On cyclic translation patch " << pp.name()
+                                << " forcing uniform separation of "
+                                << cycpp.separationVector() << endl;
+                            const_cast<vectorField&>(cpp.separation()) =
+                                pointField(1, cycpp.separationVector());
+                        }
+                        else
+                        {
+                            const_cast<vectorField&>(cpp.separation()) =
+                                pointField
+                                (
+                                    1,
+                                    pp[pp.size()/2].centre(mesh.points())
+                                  - pp[0].centre(mesh.points())
+                                );
+                        }
+                    }
+                    else
+                    {
+                        const_cast<vectorField&>(cpp.separation())
+                        .setSize(1);
+                    }
+                    Info<< "On coupled patch " << pp.name()
+                        << " forcing uniform separation of "
+                        << cpp.separation() << endl;
+                }
+                else if (!cpp.parallel())
+                {
+                    Info<< "On coupled patch " << pp.name()
+                        << " forcing uniform rotation of "
+                        << cpp.forwardT()[0] << endl;
+
+                    const_cast<tensorField&>
+                    (
+                        cpp.forwardT()
+                    ).setSize(1);
+                    const_cast<tensorField&>
+                    (
+                        cpp.reverseT()
+                    ).setSize(1);
+
+                    Info<< "On coupled patch " << pp.name()
+                        << " forcing uniform rotation of "
+                        << cpp.forwardT() << endl;
+                }
+            }
+        }
+
         Info<< "Synchronising points." << endl;
 
         pointField newPoints(mesh.points());
-        syncTools::syncPointList
+
+        syncPoints
         (
             mesh,
             newPoints,
-            nearestEqOp(),              // cop
-            point(GREAT, GREAT, GREAT), // nullValue
-            true                        // applySeparation
+            nearestEqOp(),
+            point(GREAT, GREAT, GREAT)
         );
 
         scalarField diff(mag(newPoints-mesh.points()));
diff --git a/applications/utilities/mesh/manipulation/createPatch/createPatchDict b/applications/utilities/mesh/manipulation/createPatch/createPatchDict
index 1153d0204d..5f3597f21a 100644
--- a/applications/utilities/mesh/manipulation/createPatch/createPatchDict
+++ b/applications/utilities/mesh/manipulation/createPatch/createPatchDict
@@ -1,58 +1,77 @@
-/*--------------------------------*- C++ -*----------------------------------*\
+/*---------------------------------------------------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.5                                   |
-|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|  \\    /   O peration     | Version:  1.0                                   |
+|   \\  /    A nd           | Web:      http://www.openfoam.org               |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
+
 FoamFile
 {
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    object      createPatchDict;
+    version         2.0;
+    format          ascii;
+
+    root            "";
+    case            "";
+    instance        "system";
+    local           "";
+
+    class           dictionary;
+    object          createPatcheDict;
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 // Tolerance used in matching faces. Absolute tolerance is span of
 // face times this factor.
-matchTolerance 1E-3;
+matchTolerance 1E-6;
 
 // Do a synchronisation of coupled points.
 pointSync true;
 
+
 // Patches to create.
 // If no patches does a coupled point and face synchronisation anyway.
 patches
 (
     {
         // Name of new patch
-        name leftRight0;
+        name sidePatches;
 
-        // Type of new patch
-        type cyclic;
+        // Dictionary for new patch
+        dictionary
+        {
+            type cyclic;
+            // Optional: used when matching and synchronising points.
+            //transform translational;
+            //separationVector (-2289 0 0);
+        }
 
         // How to construct: either 'patches' or 'set'
         constructFrom patches;
 
         // If constructFrom = patches : names of patches
-        patches (half0 half1);
+        //patches (periodic-1 periodic-2);
+        patches (outlet-side1 outlet-side2);
 
         // If constructFrom = set : name of faceSet
         set f0;
     }
 
-    {
-        name bottom;
-        type patch;
-
-        constructFrom set;
-
-        patches (half0 half1);
-
-        set bottomFaces;
-    }
+    //{
+    //    name bottom;
+    //    // Dictionary for new patch
+    //    dictionary
+    //    {
+    //        type patch;
+    //    }
+    //
+    //    constructFrom set;
+    //
+    //    patches (half0 half1);
+    //
+    //    set bottomFaces;
+    //}
 
 );
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C
index f4e073d4e1..73aff8a1fd 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C
@@ -263,6 +263,7 @@ void Foam::coupledPolyPatch::calcTransformTensors
         Pout<< "coupledPolyPatch::calcTransformTensors : " << name() << endl
             << "    (half)size:" << Cf.size() << nl
             << "    absTol:" << absTol << nl
+            //<< "    smallDist:" << smallDist << nl
             << "    sum(mag(nf & nr)):" << sum(mag(nf & nr)) << endl;
     }
 
@@ -316,6 +317,13 @@ void Foam::coupledPolyPatch::calcTransformTensors
             {
                 forwardT_.setSize(1);
                 reverseT_.setSize(1);
+
+                if (debug)
+                {
+                    Pout<< "    rotation " << sum(mag(forwardT_ - forwardT_[0]))
+                        << " more than local tolerance " << error
+                        << ". Assuming uniform rotation." << endl;
+                }
             }
         }
         else
@@ -384,7 +392,7 @@ void Foam::coupledPolyPatch::calcTransformTensors
 
     if (debug)
     {
-        Pout<< "    separation_:" << separation_ << nl
+        Pout<< "    separation_:" << separation_.size() << nl
             << "    forwardT size:" << forwardT_.size() << endl;
     }
 }
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C
index 31faea5669..171fb099f1 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C
@@ -127,6 +127,22 @@ void Foam::cyclicPolyPatch::calcTransforms()
             Pout<< "cyclicPolyPatch::calcTransforms : Writing half1"
                 << " faces to OBJ file " << nm1 << endl;
             writeOBJ(nm1, half1, half1.points());
+
+            OFstream str(casePath/name()+"_half0_to_half1.obj");
+            label vertI = 0;
+            Pout<< "cyclicPolyPatch::calcTransforms :"
+                << " Writing coupled face centres as lines to " << str.name()
+                << endl;
+            forAll(half0Ctrs, i)
+            {
+                const point& p0 = half0Ctrs[i];
+                str << "v " << p0.x() << ' ' << p0.y() << ' ' << p0.z() << nl;
+                vertI++;
+                const point& p1 = half1Ctrs[i];
+                str << "v " << p1.x() << ' ' << p1.y() << ' ' << p1.z() << nl;
+                vertI++;
+                str << "l " << vertI-1 << ' ' << vertI << nl;
+            }
         }
 
         vectorField half0Normals(half0.size());
@@ -397,8 +413,6 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
     anchors0 = getAnchorPoints(half0Faces, pp.points());
     half1Ctrs = calcFaceCentres(half1Faces, pp.points());
 
-    vector n0 = vector::zero;
-    vector n1 = vector::zero;
     switch (transform_)
     {
         case ROTATIONAL:
@@ -406,12 +420,46 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
             label face0 = getConsistentRotationFace(half0Ctrs);
             label face1 = getConsistentRotationFace(half1Ctrs);
 
-            n0 = ((half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_);
-            n1 = ((half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_);
+            vector n0 = ((half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_);
+            vector n1 = ((half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_);
             n0 /= mag(n0) + VSMALL;
             n1 /= mag(n1) + VSMALL;
+
+            if (debug)
+            {
+                Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
+                    << " Specified rotation :"
+                    << " n0:" << n0 << " n1:" << n1 << endl;
+            }
+
+            // Rotation (around origin)
+            const tensor reverseT(rotationTensor(n0, -n1));
+
+            // Rotation
+            forAll(half0Ctrs, faceI)
+            {
+                half0Ctrs[faceI] = Foam::transform(reverseT, half0Ctrs[faceI]);
+                anchors0[faceI] = Foam::transform(reverseT, anchors0[faceI]);
+            }
+
             break;
         }
+        //- Problem: usually specified translation is not accurate enough
+        //- to get proper match so keep automatic determination over here.
+        //case TRANSLATIONAL:
+        //{
+        //    // Transform 0 points.
+        //
+        //    if (debug)
+        //    {
+        //        Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
+        //            << "Specified translation : " << separationVector_ << endl;
+        //    }
+        //
+        //    half0Ctrs += separationVector_;
+        //    anchors0 += separationVector_;
+        //    break;
+        //}
         default:
         {
             // Assumes that cyclic is planar. This is also the initial
@@ -420,55 +468,67 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
             // Determine the face with max area on both halves. These
             // two faces are used to determine the transformation tensors
             label max0I = findMaxArea(pp.points(), half0Faces);
-            n0 = half0Faces[max0I].normal(pp.points());
+            vector n0 = half0Faces[max0I].normal(pp.points());
             n0 /= mag(n0) + VSMALL;
 
             label max1I = findMaxArea(pp.points(), half1Faces);
-            n1 = half1Faces[max1I].normal(pp.points());
+            vector n1 = half1Faces[max1I].normal(pp.points());
             n1 /= mag(n1) + VSMALL;
+
+            if (mag(n0 & n1) < 1-coupledPolyPatch::matchTol)
+            {
+                if (debug)
+                {
+                    Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
+                        << " Detected rotation :"
+                        << " n0:" << n0 << " n1:" << n1 << endl;
+                }
+
+                // Rotation (around origin)
+                const tensor reverseT(rotationTensor(n0, -n1));
+
+                // Rotation
+                forAll(half0Ctrs, faceI)
+                {
+                    half0Ctrs[faceI] = Foam::transform
+                    (
+                        reverseT,
+                        half0Ctrs[faceI]
+                    );
+                    anchors0[faceI] = Foam::transform
+                    (
+                        reverseT,
+                        anchors0[faceI]
+                    );
+                }
+            }
+            else
+            {
+                // Parallel translation. Get average of all used points.
+
+                primitiveFacePatch half0(half0Faces, pp.points());
+                const pointField& half0Pts = half0.localPoints();
+                const point ctr0(sum(half0Pts)/half0Pts.size());
+
+                primitiveFacePatch half1(half1Faces, pp.points());
+                const pointField& half1Pts = half1.localPoints();
+                const point ctr1(sum(half1Pts)/half1Pts.size());
+
+                if (debug)
+                {
+                    Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
+                        << " Detected translation :"
+                        << " n0:" << n0 << " n1:" << n1
+                        << " ctr0:" << ctr0 << " ctr1:" << ctr1 << endl;
+                }
+
+                half0Ctrs += ctr1 - ctr0;
+                anchors0 += ctr1 - ctr0;
+            }
+            break;
         }
     }
 
-    if (mag(n0 & n1) < 1-coupledPolyPatch::matchTol)
-    {
-        if (debug)
-        {
-            Pout<< "cyclicPolyPatch::getCentresAndAnchors : Rotation :"
-                << " n0:" << n0 << " n1:" << n1 << endl;
-        }
-
-        // Rotation (around origin)
-        const tensor reverseT(rotationTensor(n0, -n1));
-
-        // Rotation
-        forAll(half0Ctrs, faceI)
-        {
-            half0Ctrs[faceI] = Foam::transform(reverseT, half0Ctrs[faceI]);
-            anchors0[faceI] = Foam::transform(reverseT, anchors0[faceI]);
-        }
-    }
-    else
-    {
-        // Parallel translation. Get average of all used points.
-
-        primitiveFacePatch half0(half0Faces, pp.points());
-        const pointField& half0Pts = half0.localPoints();
-        const point ctr0(sum(half0Pts)/half0Pts.size());
-
-        primitiveFacePatch half1(half1Faces, pp.points());
-        const pointField& half1Pts = half1.localPoints();
-        const point ctr1(sum(half1Pts)/half1Pts.size());
-
-        if (debug)
-        {
-            Pout<< "cyclicPolyPatch::getCentresAndAnchors : Translation :"
-                << " n0:" << n0 << " n1:" << n1
-                << " ctr0:" << ctr0 << " ctr1:" << ctr1 << endl;
-        }
-
-        half0Ctrs += ctr1 - ctr0;
-        anchors0 += ctr1 - ctr0;
-    }
 
     // Calculate typical distance per face
     tols = calcFaceTol(half1Faces, pp.points(), half1Ctrs);
@@ -615,7 +675,8 @@ Foam::cyclicPolyPatch::cyclicPolyPatch
     featureCos_(0.9),
     transform_(UNKNOWN),
     rotationAxis_(vector::zero),
-    rotationCentre_(point::zero)
+    rotationCentre_(point::zero),
+    separationVector_(vector::zero)
 {
     calcTransforms();
 }
@@ -635,7 +696,8 @@ Foam::cyclicPolyPatch::cyclicPolyPatch
     featureCos_(0.9),
     transform_(UNKNOWN),
     rotationAxis_(vector::zero),
-    rotationCentre_(point::zero)
+    rotationCentre_(point::zero),
+    separationVector_(vector::zero)
 {
     dict.readIfPresent("featureCos", featureCos_);
 
@@ -650,9 +712,14 @@ Foam::cyclicPolyPatch::cyclicPolyPatch
                 dict.lookup("rotationCentre") >> rotationCentre_;
                 break;
             }
+            case TRANSLATIONAL:
+            {
+                dict.lookup("separationVector") >> separationVector_;
+                break;
+            }
             default:
             {
-                // no additioanl info required
+                // no additional info required
             }
         }
     }
@@ -673,7 +740,8 @@ Foam::cyclicPolyPatch::cyclicPolyPatch
     featureCos_(pp.featureCos_),
     transform_(pp.transform_),
     rotationAxis_(pp.rotationAxis_),
-    rotationCentre_(pp.rotationCentre_)
+    rotationCentre_(pp.rotationCentre_),
+    separationVector_(pp.separationVector_)
 {
     calcTransforms();
 }
@@ -694,7 +762,8 @@ Foam::cyclicPolyPatch::cyclicPolyPatch
     featureCos_(pp.featureCos_),
     transform_(pp.transform_),
     rotationAxis_(pp.rotationAxis_),
-    rotationCentre_(pp.rotationCentre_)
+    rotationCentre_(pp.rotationCentre_),
+    separationVector_(pp.separationVector_)
 {
     calcTransforms();
 }
@@ -1322,6 +1391,8 @@ void Foam::cyclicPolyPatch::write(Ostream& os) const
         {
             os.writeKeyword("transform") << transformTypeNames[TRANSLATIONAL]
                 << token::END_STATEMENT << nl;
+            os.writeKeyword("separationVector") << separationVector_
+                << token::END_STATEMENT << nl;
             break;
         }
         default:
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H
index f9c4e66dbb..7db828985f 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H
@@ -101,11 +101,18 @@ private:
         //- Type of transformation - rotational or translational
         transformType transform_;
 
-        //- Axis of rotation for rotational cyclics
-        vector rotationAxis_;
+        // For rotation
 
-        //- point on axis of rotation for rotational cyclics
-        point rotationCentre_;
+            //- Axis of rotation for rotational cyclics
+            vector rotationAxis_;
+
+            //- point on axis of rotation for rotational cyclics
+            point rotationCentre_;
+
+        // For translation
+
+            //- Translation vector
+            vector separationVector_;
 
 
     // Private member functions
@@ -267,66 +274,94 @@ public:
         const edgeList& coupledEdges() const;
 
 
-        vector separation(const label facei) const
-        {
-            if (facei < size()/2)
-            {
-                return coupledPolyPatch::separation()[0];
-            }
-            else
-            {
-                return -coupledPolyPatch::separation()[0];
-            }
-        }
 
-        const tensor& transformT(const label facei) const
-        {
-            if (facei < size()/2)
-            {
-                return reverseT()[0];
-            }
-            else
-            {
-                return forwardT()[0];
-            }
-        }
+        // Transformation
 
-        template<class T>
-        T transform(const T& t, const label facei) const
-        {
-            if (parallel())
+            vector separation(const label facei) const
             {
-                return t;
+                if (facei < size()/2)
+                {
+                    return coupledPolyPatch::separation()[0];
+                }
+                else
+                {
+                    return -coupledPolyPatch::separation()[0];
+                }
             }
-            else
-            {
-                return Foam::transform(transformT(facei), t);
-            }
-        }
 
-        label transformLocalFace(const label facei) const
-        {
-            if (facei < size()/2)
+            const tensor& transformT(const label facei) const
             {
-                return facei + size()/2;
+                if (facei < size()/2)
+                {
+                    return reverseT()[0];
+                }
+                else
+                {
+                    return forwardT()[0];
+                }
             }
-            else
-            {
-                return facei - size()/2;
-            }
-        }
 
-        label transformGlobalFace(const label facei) const
-        {
-            if (facei - start() < size()/2)
+            template<class T>
+            T transform(const T& t, const label facei) const
             {
-                return facei + size()/2;
+                if (parallel())
+                {
+                    return t;
+                }
+                else
+                {
+                    return Foam::transform(transformT(facei), t);
+                }
             }
-            else
+
+            label transformLocalFace(const label facei) const
             {
-                return facei - size()/2;
+                if (facei < size()/2)
+                {
+                    return facei + size()/2;
+                }
+                else
+                {
+                    return facei - size()/2;
+                }
             }
-        }
+
+            label transformGlobalFace(const label facei) const
+            {
+                if (facei - start() < size()/2)
+                {
+                    return facei + size()/2;
+                }
+                else
+                {
+                    return facei - size()/2;
+                }
+            }
+
+            //- Type of transform
+            transformType transform() const
+            {
+                return transform_;
+            }
+
+            //- Axis of rotation for rotational cyclics
+            const vector& rotationAxis() const
+            {
+                return rotationAxis_;
+            }
+
+            //- point on axis of rotation for rotational cyclics
+            const point& rotationCentre() const
+            {
+                return rotationCentre_;
+            }
+
+            //- Translation vector for translational cyclics
+            const vector& separationVector() const
+            {
+                return separationVector_;
+            }
+
 
 
         //- Initialize ordering for primitivePatch. Does not

From 45422cda59e62450cbf6b4d0ce3c61ba1cd2d230 Mon Sep 17 00:00:00 2001
From: henry <Henry Weller h.weller@opencfd.co.uk>
Date: Thu, 9 Oct 2008 11:23:05 +0100
Subject: [PATCH 09/10] Changed the access to derived properties to ensure the
 old-time level values are stored if required.

---
 .../basic/hThermo/hThermo.C                   | 54 ++++++------
 .../hMixtureThermo/hMixtureThermo.C           | 52 ++++++------
 .../hhuMixtureThermo/hhuMixtureThermo.C       | 82 +++++++++++--------
 3 files changed, 101 insertions(+), 87 deletions(-)

diff --git a/src/thermophysicalModels/basic/hThermo/hThermo.C b/src/thermophysicalModels/basic/hThermo/hThermo.C
index 19e43d2410..7e32f7a063 100644
--- a/src/thermophysicalModels/basic/hThermo/hThermo.C
+++ b/src/thermophysicalModels/basic/hThermo/hThermo.C
@@ -28,15 +28,10 @@ License
 #include "fvMesh.H"
 #include "fixedValueFvPatchFields.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class MixtureType>
-hThermo<MixtureType>::hThermo(const fvMesh& mesh)
+Foam::hThermo<MixtureType>::hThermo(const fvMesh& mesh)
 :
     basicThermo(mesh),
     MixtureType(*this, mesh),
@@ -56,9 +51,12 @@ hThermo<MixtureType>::hThermo(const fvMesh& mesh)
         hBoundaryTypes()
     )
 {
-    forAll(h_, celli)
+    scalarField& hCells = h_.internalField();
+    const scalarField& TCells = T_.internalField();
+
+    forAll(hCells, celli)
     {
-        h_[celli] = this->cellMixture(celli).H(T_[celli]);
+        hCells[celli] = this->cellMixture(celli).H(TCells[celli]);
     }
 
     forAll(h_.boundaryField(), patchi)
@@ -76,25 +74,33 @@ hThermo<MixtureType>::hThermo(const fvMesh& mesh)
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 template<class MixtureType>
-hThermo<MixtureType>::~hThermo()
+Foam::hThermo<MixtureType>::~hThermo()
 {}
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<class MixtureType>
-void hThermo<MixtureType>::calculate()
+void Foam::hThermo<MixtureType>::calculate()
 {
-    forAll(T_, celli)
+    const scalarField& hCells = h_.internalField();
+    const scalarField& pCells = p_.internalField();
+
+    scalarField& TCells = T_.internalField();
+    scalarField& psiCells = psi_.internalField();
+    scalarField& muCells = mu_.internalField();
+    scalarField& alphaCells = alpha_.internalField();
+
+    forAll(TCells, celli)
     {
         const typename MixtureType::thermoType& mixture_ =
             this->cellMixture(celli);
 
-        T_[celli] = mixture_.TH(h_[celli], T_[celli]);
-        psi_[celli] = mixture_.psi(p_[celli], T_[celli]);
+        TCells[celli] = mixture_.TH(hCells[celli], TCells[celli]);
+        psiCells[celli] = mixture_.psi(pCells[celli], TCells[celli]);
 
-        mu_[celli] = mixture_.mu(T_[celli]);
-        alpha_[celli] = mixture_.alpha(T_[celli]);
+        muCells[celli] = mixture_.mu(TCells[celli]);
+        alphaCells[celli] = mixture_.alpha(TCells[celli]);
     }
 
     forAll(T_.boundaryField(), patchi)
@@ -143,7 +149,7 @@ void hThermo<MixtureType>::calculate()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class MixtureType>
-void hThermo<MixtureType>::correct()
+void Foam::hThermo<MixtureType>::correct()
 {
     if (debug)
     {
@@ -163,7 +169,7 @@ void hThermo<MixtureType>::correct()
 
 
 template<class MixtureType>
-tmp<scalarField>  hThermo<MixtureType>::h
+Foam::tmp<Foam::scalarField> Foam::hThermo<MixtureType>::h
 (
     const scalarField& T,
     const labelList& cells
@@ -182,7 +188,7 @@ tmp<scalarField>  hThermo<MixtureType>::h
 
 
 template<class MixtureType>
-tmp<scalarField>  hThermo<MixtureType>::h
+Foam::tmp<Foam::scalarField> Foam::hThermo<MixtureType>::h
 (
     const scalarField& T,
     const label patchi
@@ -201,7 +207,7 @@ tmp<scalarField>  hThermo<MixtureType>::h
 
 
 template<class MixtureType>
-tmp<scalarField> hThermo<MixtureType>::Cp
+Foam::tmp<Foam::scalarField> Foam::hThermo<MixtureType>::Cp
 (
     const scalarField& T,
     const label patchi
@@ -220,7 +226,7 @@ tmp<scalarField> hThermo<MixtureType>::Cp
 
 
 template<class MixtureType>
-tmp<volScalarField> hThermo<MixtureType>::Cp() const
+Foam::tmp<Foam::volScalarField> Foam::hThermo<MixtureType>::Cp() const
 {
     const fvMesh& mesh = T_.mesh();
 
@@ -258,7 +264,7 @@ tmp<volScalarField> hThermo<MixtureType>::Cp() const
 
 
 template<class MixtureType>
-tmp<volScalarField> hThermo<MixtureType>::Cv() const
+Foam::tmp<Foam::volScalarField> Foam::hThermo<MixtureType>::Cv() const
 {
     const fvMesh& mesh = T_.mesh();
 
@@ -303,7 +309,7 @@ tmp<volScalarField> hThermo<MixtureType>::Cv() const
 
 
 template<class MixtureType>
-bool hThermo<MixtureType>::read()
+bool Foam::hThermo<MixtureType>::read()
 {
     if (basicThermo::read())
     {
@@ -317,8 +323,4 @@ bool hThermo<MixtureType>::read()
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/combustion/mixtureThermos/hMixtureThermo/hMixtureThermo.C b/src/thermophysicalModels/combustion/mixtureThermos/hMixtureThermo/hMixtureThermo.C
index 2d72a4f688..064e7c3eba 100644
--- a/src/thermophysicalModels/combustion/mixtureThermos/hMixtureThermo/hMixtureThermo.C
+++ b/src/thermophysicalModels/combustion/mixtureThermos/hMixtureThermo/hMixtureThermo.C
@@ -28,22 +28,20 @@ License
 #include "fvMesh.H"
 #include "fixedValueFvPatchFields.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class MixtureType>
-hMixtureThermo<MixtureType>::hMixtureThermo(const fvMesh& mesh)
+Foam::hMixtureThermo<MixtureType>::hMixtureThermo(const fvMesh& mesh)
 :
     hCombustionThermo(mesh),
     MixtureType(*this, mesh)
 {
-    forAll(h_, celli)
+    scalarField& hCells = h_.internalField();
+    const scalarField& TCells = T_.internalField();
+
+    forAll(hCells, celli)
     {
-        h_[celli] = this->cellMixture(celli).H(T_[celli]);
+        hCells[celli] = this->cellMixture(celli).H(TCells[celli]);
     }
 
     forAll(h_.boundaryField(), patchi)
@@ -61,25 +59,33 @@ hMixtureThermo<MixtureType>::hMixtureThermo(const fvMesh& mesh)
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 template<class MixtureType>
-hMixtureThermo<MixtureType>::~hMixtureThermo()
+Foam::hMixtureThermo<MixtureType>::~hMixtureThermo()
 {}
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<class MixtureType>
-void hMixtureThermo<MixtureType>::calculate()
+void Foam::hMixtureThermo<MixtureType>::calculate()
 {
-    forAll(T_, celli)
+    const scalarField& hCells = h_.internalField();
+    const scalarField& pCells = p_.internalField();
+
+    scalarField& TCells = T_.internalField();
+    scalarField& psiCells = psi_.internalField();
+    scalarField& muCells = mu_.internalField();
+    scalarField& alphaCells = alpha_.internalField();
+
+    forAll(TCells, celli)
     {
         const typename MixtureType::thermoType& mixture_ =
             this->cellMixture(celli);
 
-        T_[celli] = mixture_.TH(h_[celli], T_[celli]);
-        psi_[celli] = mixture_.psi(p_[celli], T_[celli]);
+        TCells[celli] = mixture_.TH(hCells[celli], TCells[celli]);
+        psiCells[celli] = mixture_.psi(pCells[celli], TCells[celli]);
 
-        mu_[celli] = mixture_.mu(T_[celli]);
-        alpha_[celli] = mixture_.alpha(T_[celli]);
+        muCells[celli] = mixture_.mu(TCells[celli]);
+        alphaCells[celli] = mixture_.alpha(TCells[celli]);
     }
 
     forAll(T_.boundaryField(), patchi)
@@ -128,7 +134,7 @@ void hMixtureThermo<MixtureType>::calculate()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class MixtureType>
-tmp<scalarField> hMixtureThermo<MixtureType>::h
+Foam::tmp<Foam::scalarField> Foam::hMixtureThermo<MixtureType>::h
 (
     const scalarField& T,
     const labelList& cells
@@ -147,7 +153,7 @@ tmp<scalarField> hMixtureThermo<MixtureType>::h
 
 
 template<class MixtureType>
-tmp<scalarField> hMixtureThermo<MixtureType>::h
+Foam::tmp<Foam::scalarField> Foam::hMixtureThermo<MixtureType>::h
 (
     const scalarField& T,
     const label patchi
@@ -166,7 +172,7 @@ tmp<scalarField> hMixtureThermo<MixtureType>::h
 
 
 template<class MixtureType>
-tmp<scalarField> hMixtureThermo<MixtureType>::Cp
+Foam::tmp<Foam::scalarField> Foam::hMixtureThermo<MixtureType>::Cp
 (
     const scalarField& T,
     const label patchi
@@ -186,7 +192,7 @@ tmp<scalarField> hMixtureThermo<MixtureType>::Cp
 
 
 template<class MixtureType>
-tmp<volScalarField> hMixtureThermo<MixtureType>::Cp() const
+Foam::tmp<Foam::volScalarField> Foam::hMixtureThermo<MixtureType>::Cp() const
 {
     const fvMesh& mesh = T_.mesh();
 
@@ -224,7 +230,7 @@ tmp<volScalarField> hMixtureThermo<MixtureType>::Cp() const
 
 
 template<class MixtureType>
-void hMixtureThermo<MixtureType>::correct()
+void Foam::hMixtureThermo<MixtureType>::correct()
 {
     if (debug)
     {
@@ -244,7 +250,7 @@ void hMixtureThermo<MixtureType>::correct()
 
 
 template<class MixtureType>
-bool hMixtureThermo<MixtureType>::read()
+bool Foam::hMixtureThermo<MixtureType>::read()
 {
     if (hCombustionThermo::read())
     {
@@ -258,8 +264,4 @@ bool hMixtureThermo<MixtureType>::read()
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/combustion/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.C b/src/thermophysicalModels/combustion/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.C
index 04bf08e330..0086e8cd14 100644
--- a/src/thermophysicalModels/combustion/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.C
+++ b/src/thermophysicalModels/combustion/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.C
@@ -28,23 +28,23 @@ License
 #include "fvMesh.H"
 #include "fixedValueFvPatchFields.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class MixtureType>
-hhuMixtureThermo<MixtureType>::hhuMixtureThermo(const fvMesh& mesh)
+Foam::hhuMixtureThermo<MixtureType>::hhuMixtureThermo(const fvMesh& mesh)
 :
     hhuCombustionThermo(mesh),
     MixtureType(*this, mesh)
 {
-    forAll(h_, celli)
+    scalarField& hCells = h_.internalField();
+    scalarField& huCells = hu_.internalField();
+    const scalarField& TCells = T_.internalField();
+    const scalarField& TuCells = Tu_.internalField();
+
+    forAll(hCells, celli)
     {
-        h_[celli] = this->cellMixture(celli).H(T_[celli]);
-        hu_[celli] = this->cellReactants(celli).H(Tu_[celli]);
+        hCells[celli] = this->cellMixture(celli).H(TCells[celli]);
+        huCells[celli] = this->cellReactants(celli).H(TuCells[celli]);
     }
 
     forAll(h_.boundaryField(), patchi)
@@ -71,27 +71,38 @@ hhuMixtureThermo<MixtureType>::hhuMixtureThermo(const fvMesh& mesh)
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 template<class MixtureType>
-hhuMixtureThermo<MixtureType>::~hhuMixtureThermo()
+Foam::hhuMixtureThermo<MixtureType>::~hhuMixtureThermo()
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class MixtureType>
-void hhuMixtureThermo<MixtureType>::calculate()
+void Foam::hhuMixtureThermo<MixtureType>::calculate()
 {
-    forAll(T_, celli)
+    const scalarField& hCells = h_.internalField();
+    const scalarField& huCells = hu_.internalField();
+    const scalarField& pCells = p_.internalField();
+
+    scalarField& TCells = T_.internalField();
+    scalarField& TuCells = Tu_.internalField();
+    scalarField& psiCells = psi_.internalField();
+    scalarField& muCells = mu_.internalField();
+    scalarField& alphaCells = alpha_.internalField();
+
+    forAll(TCells, celli)
     {
-        const typename MixtureType::thermoType& mixture_ = 
+        const typename MixtureType::thermoType& mixture_ =
             this->cellMixture(celli);
 
-        T_[celli] = mixture_.TH(h_[celli], T_[celli]);
-        psi_[celli] = mixture_.psi(p_[celli], T_[celli]);
+        TCells[celli] = mixture_.TH(hCells[celli], TCells[celli]);
+        psiCells[celli] = mixture_.psi(pCells[celli], TCells[celli]);
 
-        mu_[celli] = mixture_.mu(T_[celli]);
-        alpha_[celli] = mixture_.alpha(T_[celli]);
+        muCells[celli] = mixture_.mu(TCells[celli]);
+        alphaCells[celli] = mixture_.alpha(TCells[celli]);
 
-        Tu_[celli] = this->cellReactants(celli).TH(hu_[celli], Tu_[celli]);
+        TuCells[celli] =
+            this->cellReactants(celli).TH(huCells[celli], TuCells[celli]);
     }
 
     forAll(T_.boundaryField(), patchi)
@@ -144,7 +155,7 @@ void hhuMixtureThermo<MixtureType>::calculate()
 
 
 template<class MixtureType>
-void hhuMixtureThermo<MixtureType>::correct()
+void Foam::hhuMixtureThermo<MixtureType>::correct()
 {
     if (debug)
     {
@@ -163,7 +174,7 @@ void hhuMixtureThermo<MixtureType>::correct()
 }
 
 template<class MixtureType>
-tmp<scalarField> hhuMixtureThermo<MixtureType>::h
+Foam::tmp<Foam::scalarField> Foam::hhuMixtureThermo<MixtureType>::h
 (
     const scalarField& T,
     const labelList& cells
@@ -182,7 +193,7 @@ tmp<scalarField> hhuMixtureThermo<MixtureType>::h
 
 
 template<class MixtureType>
-tmp<scalarField> hhuMixtureThermo<MixtureType>::h
+Foam::tmp<Foam::scalarField> Foam::hhuMixtureThermo<MixtureType>::h
 (
     const scalarField& T,
     const label patchi
@@ -201,7 +212,7 @@ tmp<scalarField> hhuMixtureThermo<MixtureType>::h
 
 
 template<class MixtureType>
-tmp<scalarField> hhuMixtureThermo<MixtureType>::Cp
+Foam::tmp<Foam::scalarField> Foam::hhuMixtureThermo<MixtureType>::Cp
 (
     const scalarField& T,
     const label patchi
@@ -221,7 +232,7 @@ tmp<scalarField> hhuMixtureThermo<MixtureType>::Cp
 
 
 template<class MixtureType>
-tmp<volScalarField> hhuMixtureThermo<MixtureType>::Cp() const
+Foam::tmp<Foam::volScalarField> Foam::hhuMixtureThermo<MixtureType>::Cp() const
 {
     const fvMesh& mesh = T_.mesh();
 
@@ -260,7 +271,7 @@ tmp<volScalarField> hhuMixtureThermo<MixtureType>::Cp() const
 
 
 template<class MixtureType>
-tmp<scalarField> hhuMixtureThermo<MixtureType>::hu
+Foam::tmp<Foam::scalarField> Foam::hhuMixtureThermo<MixtureType>::hu
 (
     const scalarField& Tu,
     const labelList& cells
@@ -279,7 +290,7 @@ tmp<scalarField> hhuMixtureThermo<MixtureType>::hu
 
 
 template<class MixtureType>
-tmp<scalarField> hhuMixtureThermo<MixtureType>::hu
+Foam::tmp<Foam::scalarField> Foam::hhuMixtureThermo<MixtureType>::hu
 (
     const scalarField& Tu,
     const label patchi
@@ -298,7 +309,7 @@ tmp<scalarField> hhuMixtureThermo<MixtureType>::hu
 
 
 template<class MixtureType>
-tmp<volScalarField> hhuMixtureThermo<MixtureType>::Tb() const
+Foam::tmp<Foam::volScalarField> Foam::hhuMixtureThermo<MixtureType>::Tb() const
 {
     tmp<volScalarField> tTb
     (
@@ -342,7 +353,8 @@ tmp<volScalarField> hhuMixtureThermo<MixtureType>::Tb() const
 
 
 template<class MixtureType>
-tmp<volScalarField> hhuMixtureThermo<MixtureType>::psiu() const
+Foam::tmp<Foam::volScalarField>
+Foam::hhuMixtureThermo<MixtureType>::psiu() const
 {
     tmp<volScalarField> tpsiu
     (
@@ -388,7 +400,8 @@ tmp<volScalarField> hhuMixtureThermo<MixtureType>::psiu() const
 
 
 template<class MixtureType>
-tmp<volScalarField> hhuMixtureThermo<MixtureType>::psib() const
+Foam::tmp<Foam::volScalarField>
+Foam::hhuMixtureThermo<MixtureType>::psib() const
 {
     tmp<volScalarField> tpsib
     (
@@ -426,7 +439,8 @@ tmp<volScalarField> hhuMixtureThermo<MixtureType>::psib() const
         forAll(ppsib, facei)
         {
             ppsib[facei] =
-                this->patchFaceReactants(patchi, facei).psi(pp[facei], pTb[facei]);
+                this->patchFaceReactants
+                (patchi, facei).psi(pp[facei], pTb[facei]);
         }
     }
 
@@ -435,7 +449,7 @@ tmp<volScalarField> hhuMixtureThermo<MixtureType>::psib() const
 
 
 template<class MixtureType>
-tmp<volScalarField> hhuMixtureThermo<MixtureType>::muu() const
+Foam::tmp<Foam::volScalarField> Foam::hhuMixtureThermo<MixtureType>::muu() const
 {
     tmp<volScalarField> tmuu
     (
@@ -478,7 +492,7 @@ tmp<volScalarField> hhuMixtureThermo<MixtureType>::muu() const
 
 
 template<class MixtureType>
-tmp<volScalarField> hhuMixtureThermo<MixtureType>::mub() const
+Foam::tmp<Foam::volScalarField> Foam::hhuMixtureThermo<MixtureType>::mub() const
 {
     tmp<volScalarField> tmub
     (
@@ -521,7 +535,7 @@ tmp<volScalarField> hhuMixtureThermo<MixtureType>::mub() const
 
 
 template<class MixtureType>
-bool hhuMixtureThermo<MixtureType>::read()
+bool Foam::hhuMixtureThermo<MixtureType>::read()
 {
     if (hhuCombustionThermo::read())
     {
@@ -535,8 +549,4 @@ bool hhuMixtureThermo<MixtureType>::read()
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //

From b4263cca93beac0a5996dc6cfc7b8dba723581cd Mon Sep 17 00:00:00 2001
From: henry <Henry Weller h.weller@opencfd.co.uk>
Date: Thu, 9 Oct 2008 14:19:24 +0100
Subject: [PATCH 10/10] Added more aggressive limiting of the quadratic
 correction.

---
 .../schemes/quadraticFit/quadraticFitData.C   | 76 ++++++++++++++-----
 1 file changed, 57 insertions(+), 19 deletions(-)

diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C
index 7271bad0ff..79f9107e97 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C
@@ -113,7 +113,9 @@ Foam::quadraticFitData::quadraticFitData
     {
         interpPolySize[faci] = calcFit(stencilPoints[faci], faci);
     }
+
     interpPolySize.write();
+
     if (debug)
     {
         Info<< "quadraticFitData::quadraticFitData() :"
@@ -206,7 +208,7 @@ Foam::label Foam::quadraticFitData::calcFit
     // calculate the matrix of the polynomial components
     scalarRectangularMatrix B(C.size(), minSize_, scalar(0));
 
-        for(label ip = 0; ip < C.size(); ip++)
+    for(label ip = 0; ip < C.size(); ip++)
     {
         const point& p = C[ip];
 
@@ -228,6 +230,7 @@ Foam::label Foam::quadraticFitData::calcFit
         pz /= scale;
 
         label is = 0;
+
         B[ip][is++] = wts[0]*wts[ip];
         B[ip][is++] = wts[0]*wts[ip]*px;
         B[ip][is++] = wts[ip]*sqr(px);
@@ -253,21 +256,36 @@ Foam::label Foam::quadraticFitData::calcFit
     scalarList singVals(minSize_);
     label nSVDzeros = 0;
 
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>& w =
+        mesh().surfaceInterpolation::weights();
+
     bool goodFit = false;
-    for(int iIt = 0; iIt < 35 && !goodFit; iIt++)
+    for(int iIt = 0; iIt < 10 && !goodFit; iIt++)
     {
         SVD svd(B, SMALL);
 
         scalar fit0 = wts[0]*wts[0]*svd.VSinvUt()[0][0];
         scalar fit1 = wts[0]*wts[1]*svd.VSinvUt()[0][1];
 
-        goodFit = sign(fit0) == sign(fit1) && fit0 < 1 && fit1 < 1;
+        //goodFit = (fit0 > 0 && fit1 > 0);
+
+        goodFit =
+            (mag(fit0 - w[faci])/w[faci] < 0.5)
+         && (mag(fit1 - (1 - w[faci]))/(1 - w[faci]) < 0.5);
+
+        //scalar w0Err = fit0/w[faci];
+        //scalar w1Err = fit1/(1 - w[faci]);
+
+        //goodFit =
+        //    (w0Err > 0.5 && w0Err < 1.5)
+        // && (w1Err > 0.5 && w1Err < 1.5);
 
         if (goodFit)
         {
             fit_[faci][0] = fit0;
             fit_[faci][1] = fit1;
-            for(label i = 2; i < stencilSize; i++)
+
+            for(label i=2; i<stencilSize; i++)
             {
                 fit_[faci][i] = wts[0]*wts[i]*svd.VSinvUt()[0][i];
             }
@@ -279,11 +297,13 @@ Foam::label Foam::quadraticFitData::calcFit
         {
             wts[0] *= 10;
             wts[1] *= 10;
+
             for(label i = 0; i < B.n(); i++)
             {
                 B[i][0] *= 10;
                 B[i][1] *= 10;
             }
+
             for(label j = 0; j < B.m(); j++)
             {
                 B[0][j] *= 10;
@@ -291,23 +311,41 @@ Foam::label Foam::quadraticFitData::calcFit
             }
         }
     }
-    if (!goodFit)
+
+    //static const scalar alpha = 1.5;
+    //static const scalar beta = alpha/0.5;
+
+    if (goodFit)
     {
-        FatalErrorIn
-        (
-            "quadraticFitData::calcFit(const pointField&, const label)"
-            ) << "For face " << faci << endl
-            << "Cannot find good fit"
-            << exit(FatalError);
+        // scalar limiter =
+        // max
+        // (
+        //     min
+        //     (
+        //         min(alpha - beta*mag(fit_[faci][0] - w[faci])/w[faci], 1),
+        //         min(alpha - beta*mag(fit_[faci][1] - (1 - w[faci]))/(1 - w[faci]), 1)
+        //     ), 0
+        // );
+
+        // Remove the uncorrected linear coefficients
+        fit_[faci][0] -= w[faci];
+        fit_[faci][1] -= 1 - w[faci];
+
+        // if (limiter < 0.99)
+        // {
+        //     for(label i = 0; i < stencilSize; i++)
+        //     {
+        //         fit_[faci][i] *= limiter;
+        //     }
+        // }
+    }
+    else
+    {
+        Pout<< "Could not fit face " << faci
+            << " " << fit_[faci][0] << " " << w[faci]
+            << " " << fit_[faci][1] << " " << 1 - w[faci]<< endl;
+        fit_[faci] = 0;
     }
-
-    const GeometricField<scalar, fvsPatchField, surfaceMesh>& w =
-        mesh().surfaceInterpolation::weights();
-
-    // remove the uncorrected linear coefficients
-
-    fit_[faci][0] -= w[faci];
-    fit_[faci][1] -= 1 - w[faci];
 
     return minSize_ - nSVDzeros;
 }