From c1fb5e9021129343dea5d84d8baea375ce01632b Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@Germany>
Date: Wed, 6 Jul 2016 08:11:42 +0200
Subject: [PATCH 01/26] CONFIG: rename config file to avoid premature filtering
 of its name

---
 applications/utilities/mesh/conversion/Optional/Allwmake | 2 +-
 etc/config.sh/{libccmio => ccmio}                        | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)
 rename etc/config.sh/{libccmio => ccmio} (98%)

diff --git a/applications/utilities/mesh/conversion/Optional/Allwmake b/applications/utilities/mesh/conversion/Optional/Allwmake
index f84479c089..adf350c5c8 100755
--- a/applications/utilities/mesh/conversion/Optional/Allwmake
+++ b/applications/utilities/mesh/conversion/Optional/Allwmake
@@ -8,7 +8,7 @@ cd ${0%/*} || exit 1    # Run from this directory
 . $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
 # Get version info and arch-path
 . $WM_PROJECT_DIR/etc/config.sh/functions
-_foamSource $($WM_PROJECT_DIR/bin/foamEtcFile config.sh/libccmio)
+_foamSource $($WM_PROJECT_DIR/bin/foamEtcFile config.sh/ccmio)
 
 set -x
 
diff --git a/etc/config.sh/libccmio b/etc/config.sh/ccmio
similarity index 98%
rename from etc/config.sh/libccmio
rename to etc/config.sh/ccmio
index 1ef4205640..6bd6aa42ac 100644
--- a/etc/config.sh/libccmio
+++ b/etc/config.sh/ccmio
@@ -22,7 +22,7 @@
 #     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 #
 # File
-#     etc/config.sh/libccmio
+#     etc/config.sh/ccmio
 #
 # Description
 #     Setup file for libccmio include/libraries.

From b99515ea9597c40e7cb5f89a3081e6d7c8f9f14b Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@Germany>
Date: Thu, 7 Jul 2016 11:17:09 +0200
Subject: [PATCH 02/26] COMP: template ambiguity for 64-bit labels (issue #175)

---
 .../turbulentDFSEMInletFvPatchVectorField.C               | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
index 5b229917ac..b0df476fe0 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
@@ -80,7 +80,7 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::writeEddyOBJ() const
         //    label i2 = ((i + 1) % nPoint) + nPoint;
         //    os  << "l " << i1 << " " << i2 << nl;
         //}
-    }   
+    }
 
     {
         const Time& time = db().time();
@@ -95,7 +95,7 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::writeEddyOBJ() const
             const eddy& e = eddies_[eddyI];
             pointOffset += e.writeSurfaceOBJ(pointOffset, patchNormal_, os);
         }
-    }                 
+    }
 }
 
 
@@ -141,7 +141,7 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::writeLumleyCoeffs() const
             scalar eta = sqrt(-ii/3.0);
             os  << xi << token::TAB << eta << token::TAB
                 << ii << token::TAB << iii << endl;
-        }    
+        }
     }
 
     // After interpolation
@@ -1107,7 +1107,7 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::write(Ostream& os) const
     writeEntryIfDifferent(os, "d", 1.0, d_);
     writeEntryIfDifferent(os, "kappa", 0.41, kappa_);
     writeEntryIfDifferent(os, "perturb", 1e-5, perturb_);
-    writeEntryIfDifferent(os, "nCellPerEddy", 5, nCellPerEddy_);
+    writeEntryIfDifferent<label>(os, "nCellPerEddy", 5, nCellPerEddy_);
     writeEntryIfDifferent(os, "writeEddies", false, writeEddies_);
 
     if (!interpolateR_)

From fcb0b3a42d8067fd0efb2670e231a7cd2f5daebe Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@Germany>
Date: Thu, 7 Jul 2016 11:20:28 +0200
Subject: [PATCH 03/26] COMP: fftw needs int (not long) for its dimensionality
 (issue #175)

- explicitly use List<int> instead List<label> for API compatibility,
  even when 64-bit labels are in use.
---
 src/randomProcesses/Kmesh/Kmesh.C             |  2 +-
 src/randomProcesses/Kmesh/Kmesh.H             |  6 +++---
 src/randomProcesses/fft/fft.C                 | 10 +++++-----
 src/randomProcesses/fft/fft.H                 | 15 +++++++++------
 src/randomProcesses/noise/noiseFFT/noiseFFT.C |  2 +-
 5 files changed, 19 insertions(+), 16 deletions(-)

diff --git a/src/randomProcesses/Kmesh/Kmesh.C b/src/randomProcesses/Kmesh/Kmesh.C
index 5bc6ad1de3..362874cdb2 100644
--- a/src/randomProcesses/Kmesh/Kmesh.C
+++ b/src/randomProcesses/Kmesh/Kmesh.C
@@ -35,7 +35,7 @@ inline Foam::label Foam::Kmesh::index
     const label i,
     const label j,
     const label k,
-    const labelList& nn
+    const UList<int>& nn
 )
 {
     return (k + j*nn[2] + i*nn[1]*nn[2]);
diff --git a/src/randomProcesses/Kmesh/Kmesh.H b/src/randomProcesses/Kmesh/Kmesh.H
index d09a5894cb..417f1feacd 100644
--- a/src/randomProcesses/Kmesh/Kmesh.H
+++ b/src/randomProcesses/Kmesh/Kmesh.H
@@ -57,7 +57,7 @@ class Kmesh
         vector l_;
 
         //- Multi-dimensional addressing array
-        labelList nn_;
+        List<int> nn_;
 
         //- Maximum wavenumber
         scalar kmax_;
@@ -71,7 +71,7 @@ class Kmesh
             const label i,
             const label j,
             const label k,
-            const labelList& nn
+            const UList<int>& nn
         );
 
 
@@ -92,7 +92,7 @@ public:
             return l_;
         }
 
-        const labelList& nn() const
+        const List<int>& nn() const
         {
             return nn_;
         }
diff --git a/src/randomProcesses/fft/fft.C b/src/randomProcesses/fft/fft.C
index e8194da8dd..fd44c0109d 100644
--- a/src/randomProcesses/fft/fft.C
+++ b/src/randomProcesses/fft/fft.C
@@ -37,7 +37,7 @@ namespace Foam
 void fft::transform
 (
     complexField& field,
-    const labelList& nn,
+    const UList<int>& nn,
     transformDirection dir
 )
 {
@@ -112,7 +112,7 @@ void fft::transform
 tmp<complexField> fft::forwardTransform
 (
     const tmp<complexField>& tfield,
-    const labelList& nn
+    const UList<int>& nn
 )
 {
     tmp<complexField> tfftField(new complexField(tfield));
@@ -128,7 +128,7 @@ tmp<complexField> fft::forwardTransform
 tmp<complexField> fft::reverseTransform
 (
     const tmp<complexField>& tfield,
-    const labelList& nn
+    const UList<int>& nn
 )
 {
     tmp<complexField> tifftField(new complexField(tfield));
@@ -144,7 +144,7 @@ tmp<complexField> fft::reverseTransform
 tmp<complexVectorField> fft::forwardTransform
 (
     const tmp<complexVectorField>& tfield,
-    const labelList& nn
+    const UList<int>& nn
 )
 {
     tmp<complexVectorField> tfftVectorField
@@ -173,7 +173,7 @@ tmp<complexVectorField> fft::forwardTransform
 tmp<complexVectorField> fft::reverseTransform
 (
     const tmp<complexVectorField>& tfield,
-    const labelList& nn
+    const UList<int>& nn
 )
 {
     tmp<complexVectorField> tifftVectorField
diff --git a/src/randomProcesses/fft/fft.H b/src/randomProcesses/fft/fft.H
index dc5a0c2c69..11cc07575e 100644
--- a/src/randomProcesses/fft/fft.H
+++ b/src/randomProcesses/fft/fft.H
@@ -32,6 +32,9 @@ Description
     reverse).  The dimensionality and organisation of the array of values
     in space is supplied in the nn indexing array.
 
+Note
+    The fftw library uses int only (no longs) for its dimensionality.
+
 SourceFiles
     fft.C
 
@@ -41,7 +44,7 @@ SourceFiles
 #define fft_H
 
 #include "complexFields.H"
-#include "labelList.H"
+#include "UList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -63,7 +66,7 @@ public:
     static void transform
     (
         complexField& field,
-        const labelList& nn,
+        const UList<int>& nn,
         transformDirection fftDirection
     );
 
@@ -71,28 +74,28 @@ public:
     static tmp<complexField> forwardTransform
     (
         const tmp<complexField>& field,
-        const labelList& nn
+        const UList<int>& nn
     );
 
 
     static tmp<complexField> reverseTransform
     (
         const tmp<complexField>& field,
-        const labelList& nn
+        const UList<int>& nn
     );
 
 
     static tmp<complexVectorField> forwardTransform
     (
         const tmp<complexVectorField>& field,
-        const labelList& nn
+        const UList<int>& nn
     );
 
 
     static tmp<complexVectorField> reverseTransform
     (
         const tmp<complexVectorField>& field,
-        const labelList& nn
+        const UList<int>& nn
     );
 };
 
diff --git a/src/randomProcesses/noise/noiseFFT/noiseFFT.C b/src/randomProcesses/noise/noiseFFT/noiseFFT.C
index abaf199514..3a61c9b0dd 100644
--- a/src/randomProcesses/noise/noiseFFT/noiseFFT.C
+++ b/src/randomProcesses/noise/noiseFFT/noiseFFT.C
@@ -246,7 +246,7 @@ Foam::tmp<Foam::scalarField> Foam::noiseFFT::Pf
             fft::reverseTransform
             (
                 ReComplexField(tpn),
-                labelList(1, tpn().size())
+                List<int>(1, tpn().size())
             )
         )
     );

From 5d594d3f24cfb14fec2454340fa23658edf9f282 Mon Sep 17 00:00:00 2001
From: sergio <sergio>
Date: Thu, 7 Jul 2016 10:19:56 -0700
Subject: [PATCH 04/26] ENH: Adding support to moving mesh to
 activePressureForceBaffle

---
 ...PressureForceBaffleVelocityFvPatchVectorField.C | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.C
index 5371eda475..ad6ebca09f 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.C
@@ -332,6 +332,20 @@ void Foam::activePressureForceBaffleVelocityFvPatchVectorField::updateCoeffs()
             areaFraction = 1 - openFraction_;
         }
 
+        if (patch().boundaryMesh().mesh().moving())
+        {
+            // All areas have been recalculated already so are consistent
+            // with the new points. Note we already only do this routine
+            // once per timestep. The alternative is to use the upToDate
+            // mechanism for regIOobject + polyMesh::upToDatePoints
+            initWallSf_ = patch().Sf();
+            initCyclicSf_ = patch().boundaryMesh()[cyclicPatchLabel_].Sf();
+            nbrCyclicSf_ =  refCast<const cyclicFvPatch>
+            (
+                patch().boundaryMesh()[cyclicPatchLabel_]
+            ).neighbFvPatch().Sf();
+        }
+
         // Update this wall patch
         vectorField::subField Sfw = patch().patch().faceAreas();
         vectorField newSfw((1 - areaFraction)*initWallSf_);

From 4dcfd245d6940d43867638943c78a7ac47ea1a9f Mon Sep 17 00:00:00 2001
From: sergio <sergio>
Date: Tue, 19 Jul 2016 10:51:50 -0700
Subject: [PATCH 05/26] BUG: Taken out mesh flux from the hEqn for Solid to Gas
 reactions

---
 .../pyrolysisModels/reactingOneDim/reactingOneDim.C         | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
index af8bf5563a..2d3535d873 100644
--- a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
+++ b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
@@ -335,6 +335,10 @@ void reactingOneDim::solveEnergy()
         hEqn += fvc::div(phiQr);
     }
 
+/* NOTE: The moving mesh option is only correct for reaction such as
+   Solid -> Gas, thus the ddt term is compesated exaclty by chemistrySh and
+   the mesh flux is not necessary.
+   
     if (regionMesh().moving())
     {
         surfaceScalarField phihMesh
@@ -344,7 +348,7 @@ void reactingOneDim::solveEnergy()
 
         hEqn -= fvc::div(phihMesh);
     }
-
+*/
     hEqn.relax();
     hEqn.solve();
 }

From cd8be89b287762a06f00baf951a4096bea73cfeb Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Fri, 22 Jul 2016 15:32:17 +0100
Subject: [PATCH 06/26] BUG: runTimePostProcessing - correct the glyph
 behaviour when dealing with cell data.  Fixes #186

---
 .../fieldVisualisationBase.C                  | 140 ++++++++++--------
 1 file changed, 77 insertions(+), 63 deletions(-)

diff --git a/src/postProcessing/functionObjects/graphics/runTimePostProcessing/fieldVisualisationBase.C b/src/postProcessing/functionObjects/graphics/runTimePostProcessing/fieldVisualisationBase.C
index 40f60a7c86..9d163da19d 100644
--- a/src/postProcessing/functionObjects/graphics/runTimePostProcessing/fieldVisualisationBase.C
+++ b/src/postProcessing/functionObjects/graphics/runTimePostProcessing/fieldVisualisationBase.C
@@ -43,6 +43,7 @@ License
 #include "vtkSphereSource.h"
 #include "vtkTextActor.h"
 #include "vtkTextProperty.h"
+#include "vtkCellDataToPointData.h"
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
@@ -151,6 +152,9 @@ void Foam::fieldVisualisationBase::addScalarBar
     const vector textColour = colours_["text"]->value(position);
 
     // Work-around to supply our own scalarbar title
+    // - Default scalar bar title text is scales by the scalar bar box
+    //   dimensions so if the title is a long string, the text is shrunk to fit
+    //   Instead, suppress title and set the title using a vtkTextActor
     vtkSmartPointer<vtkTextActor> titleActor =
         vtkSmartPointer<vtkTextActor>::New();
     sbar->SetTitle(" ");
@@ -170,19 +174,18 @@ void Foam::fieldVisualisationBase::addScalarBar
     titleActor->GetPositionCoordinate()->
         SetCoordinateSystemToNormalizedViewport();
 
-/*
-    sbar->SetTitle(scalarBar_.title_.c_str());
-    sbar->GetTitleTextProperty()->SetColor
-    (
-        textColour[0],
-        textColour[1],
-        textColour[2]
-    );
-    sbar->GetTitleTextProperty()->SetFontSize(scalarBar_.fontSize_);
-    sbar->GetTitleTextProperty()->ShadowOff();
-    sbar->GetTitleTextProperty()->BoldOn();
-    sbar->GetTitleTextProperty()->ItalicOff();
-*/
+    // How to use the standard scalar bar text
+    // sbar->SetTitle(scalarBar_.title_.c_str());
+    // sbar->GetTitleTextProperty()->SetColor
+    // (
+    //     textColour[0],
+    //     textColour[1],
+    //     textColour[2]
+    // );
+    // sbar->GetTitleTextProperty()->SetFontSize(scalarBar_.fontSize_);
+    // sbar->GetTitleTextProperty()->ShadowOff();
+    // sbar->GetTitleTextProperty()->BoldOn();
+    // sbar->GetTitleTextProperty()->ItalicOff();
 
     sbar->GetLabelTextProperty()->SetColor
     (
@@ -217,8 +220,8 @@ void Foam::fieldVisualisationBase::addScalarBar
         sbar->SetWidth(0.75);
         sbar->SetHeight(0.07);
         sbar->SetBarRatio(0.5);
-//        sbar->SetHeight(0.1);
-//        sbar->SetTitleRatio(0.01);
+        // sbar->SetHeight(0.1);
+        // sbar->SetTitleRatio(0.01);
         sbar->SetTextPositionToPrecedeScalarBar();
     }
 
@@ -228,10 +231,10 @@ void Foam::fieldVisualisationBase::addScalarBar
         scalarBar_.position_.second() + sbar->GetHeight()
     );
 
-//    sbar->DrawFrameOn();
-//    sbar->DrawBackgroundOn();
-//    sbar->UseOpacityOff();
-//    sbar->VisibilityOff();
+    // sbar->DrawFrameOn();
+    // sbar->DrawBackgroundOn();
+    // sbar->UseOpacityOff();
+    // sbar->VisibilityOff();
     sbar->VisibilityOn();
 
     renderer->AddActor(sbar);
@@ -268,25 +271,7 @@ void Foam::fieldVisualisationBase::setField
             // Configure the mapper
             mapper->SelectColorArray(colourFieldName.c_str());
             mapper->SetScalarRange(range_.first(), range_.second());
-
-            // Set to use either cell or point data
-            const char* fieldName = colourFieldName.c_str();
-            if (pData->GetCellData()->HasArray(fieldName) == 1)
-            {
-                mapper->SetScalarModeToUseCellFieldData();
-            }
-            else if (pData->GetPointData()->HasArray(fieldName) == 1)
-            {
-                mapper->SetScalarModeToUsePointFieldData();
-            }
-            else
-            {
-                WarningInFunction
-                    << "Unable to determine cell or point data type "
-                    << "- assuming point data";
-                mapper->SetScalarModeToUsePointFieldData();
-            }
-
+            mapper->SetScalarModeToDefault(); // try points, then cells
             mapper->SetColorModeToMapScalars();
             mapper->SetLookupTable(lut);
             mapper->ScalarVisibilityOn();
@@ -322,9 +307,37 @@ void Foam::fieldVisualisationBase::addGlyphs
     glyph->ScalingOn();
     bool ok = true;
 
-    label nComponents =
-        data->GetPointData()->GetArray(scaleFieldName.c_str())
-            ->GetNumberOfComponents();
+    // Determine whether we have scalar or vector data
+    label nComponents = -1;
+    const char* scaleFieldNameChar = scaleFieldName.c_str();
+    if (data->GetPointData()->HasArray(scaleFieldNameChar) == 1)
+    {
+        nComponents =
+            data->GetPointData()->GetArray(scaleFieldNameChar)
+                ->GetNumberOfComponents();
+    }
+    else if (data->GetCellData()->HasArray(scaleFieldNameChar) == 1)
+    {
+        // Need to convert cell data to point data
+        vtkSmartPointer<vtkCellDataToPointData> cellToPoint =
+            vtkSmartPointer<vtkCellDataToPointData>::New();
+        cellToPoint->SetInputData(data);
+        cellToPoint->Update();
+        vtkDataSet* pds = cellToPoint->GetOutput();
+        vtkDataArray* pData = pds->GetPointData()->GetArray(scaleFieldNameChar);
+
+        // Store in main vtkPolyData
+        data->GetPointData()->AddArray(pData);
+
+        nComponents = pData->GetNumberOfComponents();
+    }
+    else
+    {
+        WarningInFunction
+            << "Glyphs can only be added to scalar or vector data. "
+            << "Unable to process field " << scaleFieldName << endl;
+        return;
+    }
 
     if (nComponents == 1)
     {
@@ -332,9 +345,10 @@ void Foam::fieldVisualisationBase::addGlyphs
             vtkSmartPointer<vtkSphereSource>::New();
         sphere->SetCenter(0, 0, 0);
         sphere->SetRadius(0.5);
-// Setting higher resolution slows the rendering significantly
-//        sphere->SetPhiResolution(20);
-//        sphere->SetThetaResolution(20);
+
+        // Setting higher resolution slows the rendering significantly
+        // sphere->SetPhiResolution(20);
+        // sphere->SetThetaResolution(20);
 
         glyph->SetSourceConnection(sphere->GetOutputPort());
 
@@ -342,18 +356,18 @@ void Foam::fieldVisualisationBase::addGlyphs
         {
             double range[2];
 
-// Can use values to find range
-//            vtkDataArray* values =
-//                data->GetPointData()->GetScalars(scaleFieldName.c_str());
-//            values->GetRange(range);
+            // Can use values to find range
+            // vtkDataArray* values =
+            //     data->GetPointData()->GetScalars(scaleFieldNameChar);
+            // values->GetRange(range);
 
-            // set range accoding to user-supplied limits
+            // Set range accoding to user-supplied limits
             range[0] = range_.first();
             range[1] = range_.second();
             glyph->ClampingOn();
             glyph->SetRange(range);
 
-            // if range[0] != min(value), maxGlyphLength behaviour will not
+            // If range[0] != min(value), maxGlyphLength behaviour will not
             // be correct...
             glyph->SetScaleFactor(maxGlyphLength);
         }
@@ -370,7 +384,7 @@ void Foam::fieldVisualisationBase::addGlyphs
             0,
             0,
             vtkDataObject::FIELD_ASSOCIATION_POINTS,
-            scaleFieldName.c_str()
+            scaleFieldNameChar
         );
     }
     else if (nComponents == 3)
@@ -388,21 +402,21 @@ void Foam::fieldVisualisationBase::addGlyphs
         if (maxGlyphLength > 0)
         {
             vtkDataArray* values =
-                data->GetPointData()->GetVectors(scaleFieldName.c_str());
+                data->GetPointData()->GetVectors(scaleFieldNameChar);
+
             double range[6];
             values->GetRange(range);
 
-/*
             // Attempt to set range for vectors...
-            scalar x0 = sqrt(sqr(range_.first())/3.0);
-            scalar x1 = sqrt(sqr(range_.second())/3.0);
-            range[0] = x0;
-            range[1] = x0;
-            range[2] = x0;
-            range[3] = x1;
-            range[4] = x1;
-            range[5] = x1;
-*/
+            // scalar x0 = sqrt(sqr(range_.first())/3.0);
+            // scalar x1 = sqrt(sqr(range_.second())/3.0);
+            // range[0] = x0;
+            // range[1] = x0;
+            // range[2] = x0;
+            // range[3] = x1;
+            // range[4] = x1;
+            // range[5] = x1;
+
             glyph->ClampingOn();
             glyph->SetRange(range);
             glyph->SetScaleFactor(maxGlyphLength);
@@ -421,7 +435,7 @@ void Foam::fieldVisualisationBase::addGlyphs
             0,
             0,
             vtkDataObject::FIELD_ASSOCIATION_POINTS,
-            scaleFieldName.c_str()
+            scaleFieldNameChar
         );
     }
     else

From 2d96fc5011de6cc5027aef78f4de646d3ca48daa Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Tue, 26 Jul 2016 12:05:15 +0100
Subject: [PATCH 07/26] BUG: mapFields function object - evaluate constraint
 patches after mapping.  Fixes #190

---
 .../utilities/mapFields/mapFields.H           |  7 ++
 .../utilities/mapFields/mapFieldsTemplates.C  | 88 +++++++++++++++++++
 2 files changed, 95 insertions(+)

diff --git a/src/postProcessing/functionObjects/utilities/mapFields/mapFields.H b/src/postProcessing/functionObjects/utilities/mapFields/mapFields.H
index 80c2ebe119..9343616cfe 100644
--- a/src/postProcessing/functionObjects/utilities/mapFields/mapFields.H
+++ b/src/postProcessing/functionObjects/utilities/mapFields/mapFields.H
@@ -122,6 +122,13 @@ class mapFieldsFO
         //- Helper function to create the mesh-to-mesh interpolation
         void createInterpolation(const dictionary& dict);
 
+        //- Helper function to evaluate constraint patches after mapping
+        template<class Type>
+        void evaluateConstraintTypes
+        (
+            GeometricField<Type, fvPatchField, volMesh>& fld
+        ) const;
+
         //- Helper function to interpolate and write the fied
         template<class Type>
         bool writeFieldType() const;
diff --git a/src/postProcessing/functionObjects/utilities/mapFields/mapFieldsTemplates.C b/src/postProcessing/functionObjects/utilities/mapFields/mapFieldsTemplates.C
index de269ea867..dd1086049b 100644
--- a/src/postProcessing/functionObjects/utilities/mapFields/mapFieldsTemplates.C
+++ b/src/postProcessing/functionObjects/utilities/mapFields/mapFieldsTemplates.C
@@ -25,6 +25,91 @@ License
 
 #include "meshToMesh.H"
 
+template<class Type>
+void Foam::mapFieldsFO::evaluateConstraintTypes
+(
+    GeometricField<Type, fvPatchField, volMesh>& fld
+) const
+{
+    typename GeometricField<Type, fvPatchField, volMesh>::
+        GeometricBoundaryField& fldBf = fld.boundaryField();
+
+    if
+    (
+        Pstream::defaultCommsType == Pstream::blocking
+     || Pstream::defaultCommsType == Pstream::nonBlocking
+    )
+    {
+        label nReq = Pstream::nRequests();
+
+        forAll(fldBf, patchi)
+        {
+            fvPatchField<Type>& tgtField = fldBf[patchi];
+
+            if
+            (
+                tgtField.type() == tgtField.patch().patch().type()
+             && polyPatch::constraintType(tgtField.patch().patch().type())
+            )
+            {
+                tgtField.initEvaluate(Pstream::defaultCommsType);
+            }
+        }
+
+        // Block for any outstanding requests
+        if
+        (
+            Pstream::parRun()
+         && Pstream::defaultCommsType == Pstream::nonBlocking
+        )
+        {
+            Pstream::waitRequests(nReq);
+        }
+
+        forAll(fldBf, patchi)
+        {
+            fvPatchField<Type>& tgtField = fldBf[patchi];
+
+            if
+            (
+                tgtField.type() == tgtField.patch().patch().type()
+             && polyPatch::constraintType(tgtField.patch().patch().type())
+            )
+            {
+                tgtField.evaluate(Pstream::defaultCommsType);
+            }
+        }
+    }
+    else if (Pstream::defaultCommsType == Pstream::scheduled)
+    {
+        const lduSchedule& patchSchedule =
+            fld.mesh().globalData().patchSchedule();
+
+        forAll(patchSchedule, patchEvali)
+        {
+            label patchi = patchSchedule[patchEvali].patch;
+            fvPatchField<Type>& tgtField = fldBf[patchi];
+
+            if
+            (
+                tgtField.type() == tgtField.patch().patch().type()
+             && polyPatch::constraintType(tgtField.patch().patch().type())
+            )
+            {
+                if (patchSchedule[patchEvali].init)
+                {
+                    tgtField.initEvaluate(Pstream::scheduled);
+                }
+                else
+                {
+                    tgtField.evaluate(Pstream::scheduled);
+                }
+            }
+        }
+    }
+}
+
+
 template<class Type>
 bool Foam::mapFieldsFO::writeFieldType() const
 {
@@ -53,6 +138,9 @@ bool Foam::mapFieldsFO::writeFieldType() const
         if (log_) Info<< ": interpolated";
 
         FieldType fieldMapRegion(mapRegionIO, tfieldMapRegion);
+
+        evaluateConstraintTypes(fieldMapRegion);
+
         fieldMapRegion.write();
 
         if (log_) Info<< " and written" << nl;

From 143800ab3d27b8d9379d8724267b4e17072e2a33 Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Tue, 26 Jul 2016 15:33:10 +0100
Subject: [PATCH 08/26] BUG: Turbulence DFSEM - protect sigFpe for reconstruct
 - see #194

---
 .../turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
index b0df476fe0..59646f3ad8 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
@@ -848,7 +848,7 @@ turbulentDFSEMInletFvPatchVectorField
     eddy::debug = debug;
 
     // Set UMean as patch area average value
-    UMean_ = gSum(U_*patch().magSf())/gSum(patch().magSf());
+    UMean_ = gSum(U_*patch().magSf())/(gSum(patch().magSf()) + ROOTVSMALL);
 }
 
 

From 785b32830d46d53e4a16af05014c3eb3b5d8becb Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Wed, 27 Jul 2016 13:06:22 +0100
Subject: [PATCH 09/26] STYLE: fanFvPatchField - updated example usage to
 reflect change in CSV specification.  Fixes #177

---
 .../fields/fvPatchFields/derived/fan/fanFvPatchField.H       | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H
index b1080e264f..8ddbaf0e32 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H
@@ -53,10 +53,11 @@ Description
         jumpTable       csvFile;
         csvFileCoeffs
         {
-            hasHeaderLine   1;
+            nHeaderLine     1;
             refColumn       0;
             componentColumns 1(1);
             separator       ",";
+            mergeSeparators no;
             fileName        "$FOAM_CASE/constant/pressureVsU";
         }
         value           uniform 0;
@@ -67,7 +68,7 @@ Description
     the jump condition.
 
 Note
-     The underlying \c patchType should be set to \c cyclic
+    The underlying \c patchType should be set to \c cyclic
 
 SeeAlso
     Foam::Function1Types

From 14f5239e5b71205378e5bea493c7bdb61c4f5c89 Mon Sep 17 00:00:00 2001
From: sergio <sergio>
Date: Wed, 27 Jul 2016 08:55:52 -0700
Subject: [PATCH 10/26] ENH: Adding optional uniform jump switch to fan BC

---
 .../derived/fan/fanFvPatchField.C             | 19 ++++++++++++++-----
 .../derived/fan/fanFvPatchField.H             |  5 +++++
 .../derived/fan/fanFvPatchFields.C            |  7 ++++++-
 3 files changed, 25 insertions(+), 6 deletions(-)

diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.C
index 3f855d52fa..24714aaf46 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.C
@@ -48,7 +48,8 @@ Foam::fanFvPatchField<Type>::fanFvPatchField
 :
     uniformJumpFvPatchField<Type>(p, iF),
     phiName_("phi"),
-    rhoName_("rho")
+    rhoName_("rho"),
+    uniformJump_(false)
 {}
 
 
@@ -62,7 +63,8 @@ Foam::fanFvPatchField<Type>::fanFvPatchField
 :
     uniformJumpFvPatchField<Type>(p, iF, dict),
     phiName_(dict.lookupOrDefault<word>("phi", "phi")),
-    rhoName_(dict.lookupOrDefault<word>("rho", "rho"))
+    rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
+    uniformJump_(dict.lookupOrDefault<bool>("uniformJump", "false"))
 {}
 
 
@@ -77,7 +79,8 @@ Foam::fanFvPatchField<Type>::fanFvPatchField
 :
     uniformJumpFvPatchField<Type>(ptf, p, iF, mapper),
     phiName_(ptf.phiName_),
-    rhoName_(ptf.rhoName_)
+    rhoName_(ptf.rhoName_),
+    uniformJump_(ptf.uniformJump_)
 {}
 
 
@@ -89,7 +92,8 @@ Foam::fanFvPatchField<Type>::fanFvPatchField
 :
     uniformJumpFvPatchField<Type>(ptf),
     phiName_(ptf.phiName_),
-    rhoName_(ptf.rhoName_)
+    rhoName_(ptf.rhoName_),
+    uniformJump_(ptf.uniformJump_)
 {}
 
 
@@ -102,7 +106,8 @@ Foam::fanFvPatchField<Type>::fanFvPatchField
 :
     uniformJumpFvPatchField<Type>(ptf, iF),
     phiName_(ptf.phiName_),
-    rhoName_(ptf.rhoName_)
+    rhoName_(ptf.rhoName_),
+    uniformJump_(ptf.uniformJump_)
 {}
 
 
@@ -129,6 +134,10 @@ void Foam::fanFvPatchField<Type>::write(Ostream& os) const
     uniformJumpFvPatchField<Type>::write(os);
     this->template writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
     this->template writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
+    this->template writeEntryIfDifferent<bool>
+    (
+        os, "uniformJump", false, uniformJump_
+    );
 }
 
 
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H
index b1080e264f..874e071159 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H
@@ -42,6 +42,8 @@ Description
         jumpTable    | jump data, e.g. \c csvFile | yes      |
         phi          | flux field name         | no          | phi
         rho          | density field name      | no          | none
+        uniformJump  | applies a uniform pressure based on the averaged
+                       velocity | no | false
     \endtable
 
     Example of the boundary condition specification:
@@ -109,6 +111,9 @@ class fanFvPatchField
         //  if neccessary
         word rhoName_;
 
+        //- Uniform pressure drop
+        bool uniformJump_;
+
 
     // Private Member Functions
 
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchFields.C b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchFields.C
index a4f2809df8..d2580b1b82 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchFields.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchFields.C
@@ -44,6 +44,10 @@ void Foam::fanFvPatchField<Foam::scalar>::calcFanJump()
             patch().patchField<surfaceScalarField, scalar>(phi);
 
         scalarField Un(max(phip/patch().magSf(), scalar(0)));
+        if (uniformJump_)
+        {
+             Un = gAverage(Un);
+        }
 
         if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
         {
@@ -67,7 +71,8 @@ Foam::fanFvPatchField<Foam::scalar>::fanFvPatchField
 :
     uniformJumpFvPatchField<scalar>(p, iF),
     phiName_(dict.lookupOrDefault<word>("phi", "phi")),
-    rhoName_(dict.lookupOrDefault<word>("rho", "rho"))
+    rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
+    uniformJump_(dict.lookupOrDefault<bool>("uniformJump", false))
 {
     if (this->cyclicPatch().owner())
     {

From e9edf89a5d3a9f9ff076bc417a4313afb95fd201 Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Mon, 1 Aug 2016 12:59:24 +0100
Subject: [PATCH 11/26] BUG: maxDeltaxyz LES delta - updated calculation. 
 Fixes #200

---
 .../LES/LESdeltas/maxDeltaxyz/maxDeltaxyz.C   | 23 +++++++++++--------
 1 file changed, 14 insertions(+), 9 deletions(-)

diff --git a/src/TurbulenceModels/turbulenceModels/LES/LESdeltas/maxDeltaxyz/maxDeltaxyz.C b/src/TurbulenceModels/turbulenceModels/LES/LESdeltas/maxDeltaxyz/maxDeltaxyz.C
index d43f4c2578..77c2d8f310 100644
--- a/src/TurbulenceModels/turbulenceModels/LES/LESdeltas/maxDeltaxyz/maxDeltaxyz.C
+++ b/src/TurbulenceModels/turbulenceModels/LES/LESdeltas/maxDeltaxyz/maxDeltaxyz.C
@@ -47,26 +47,31 @@ void Foam::LESModels::maxDeltaxyz::calcDelta()
     label nD = mesh.nGeometricD();
 
     const cellList& cells = mesh.cells();
+    const vectorField& cellC = mesh.cellCentres();
+    const vectorField& faceC = mesh.faceCentres();
+    const vectorField faceN(mesh.faceAreas()/mag(mesh.faceAreas()));
     scalarField hmax(cells.size());
 
-    forAll(cells,cellI)
+    forAll(cells, celli)
     {
         scalar deltaMaxTmp = 0.0;
-        const labelList& cFaces = mesh.cells()[cellI];
-        const point& centrevector = mesh.cellCentres()[cellI];
+        const labelList& cFaces = cells[celli];
+        const point& cc = cellC[celli];
 
-        forAll(cFaces, cFaceI)
+        forAll(cFaces, cFacei)
         {
-            label faceI = cFaces[cFaceI];
-            const point& facevector = mesh.faceCentres()[faceI];
-            scalar tmp = mag(facevector - centrevector);
+            label facei = cFaces[cFacei];
+            const point& fc = faceC[facei];
+            const vector& n = faceN[facei];
+
+            scalar tmp = magSqr(n*(n & (fc - cc)));
             if (tmp > deltaMaxTmp)
             {
                 deltaMaxTmp = tmp;
             }
         }
 
-        hmax[cellI] = deltaCoeff_*deltaMaxTmp;
+        hmax[celli] = deltaCoeff_*Foam::sqrt(deltaMaxTmp);
     }
 
     if (nD == 3)
@@ -76,7 +81,7 @@ void Foam::LESModels::maxDeltaxyz::calcDelta()
     else if (nD == 2)
     {
         WarningInFunction
-            << "Case is 2D, LES is not strictly applicable\n"
+            << "Case is 2D, LES is not strictly applicable" << nl
             << endl;
 
         delta_.internalField() = hmax;

From a29ff3f6f29bfaf64ad8edfc55f028f3044da439 Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Tue, 9 Aug 2016 15:33:46 +0100
Subject: [PATCH 12/26] ENH: boundedConvection scheme - added access function
 to the underlying scheme

---
 .../boundedConvectionScheme/boundedConvectionScheme.C      | 7 +++++++
 .../boundedConvectionScheme/boundedConvectionScheme.H      | 2 ++
 2 files changed, 9 insertions(+)

diff --git a/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionScheme.C b/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionScheme.C
index 5da9b1486c..904d317eda 100644
--- a/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionScheme.C
+++ b/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionScheme.C
@@ -40,6 +40,13 @@ namespace fv
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+template<class Type>
+const convectionScheme<Type>& boundedConvectionScheme<Type>::scheme() const
+{
+    return scheme_();
+}
+
+
 template<class Type>
 tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
 boundedConvectionScheme<Type>::interpolate
diff --git a/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionScheme.H b/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionScheme.H
index 630e17090a..0573ef4fc2 100644
--- a/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionScheme.H
+++ b/src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionScheme.H
@@ -106,6 +106,8 @@ public:
 
     // Member Functions
 
+        const convectionScheme<Type>& scheme() const;
+
         tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> interpolate
         (
             const surfaceScalarField&,

From 8509fc3ce4777912b08d5dca05f2ef9b1aaac22d Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Tue, 9 Aug 2016 15:38:40 +0100
Subject: [PATCH 13/26] BUG: blendingFactor function object updated to handle
 bounded schemes.  Fixes #205

---
 .../utilities/blendingFactor/blendingFactor.H | 12 ++-
 .../blendingFactor/blendingFactorTemplates.C  | 87 ++++++++++++++-----
 2 files changed, 74 insertions(+), 25 deletions(-)

diff --git a/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactor.H b/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactor.H
index f7166470af..be2fe747ad 100644
--- a/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactor.H
+++ b/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactor.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -88,6 +88,8 @@ SourceFiles
 #include "OFstream.H"
 #include "Switch.H"
 
+#include "convectionScheme.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -140,6 +142,14 @@ class blendingFactor
         //- Disallow default bitwise assignment
         void operator=(const blendingFactor&);
 
+        //- Helper function to calculate the blending factor for the scheme
+        template<class Type>
+        void calcScheme
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& field,
+            const typename fv::convectionScheme<Type>& cs
+        );
+
         //- Calculate the blending factor
         template<class Type>
         void calc();
diff --git a/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactorTemplates.C b/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactorTemplates.C
index 26832a892d..7031607bd4 100644
--- a/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactorTemplates.C
+++ b/src/postProcessing/functionObjects/utilities/blendingFactor/blendingFactorTemplates.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "gaussConvectionScheme.H"
+#include "boundedConvectionScheme.H"
 #include "blendedSchemeBase.H"
 #include "fvcCellReduce.H"
 #include "zeroGradientFvPatchFields.H"
@@ -31,39 +32,35 @@ License
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 template<class Type>
-void Foam::blendingFactor::calc()
+void Foam::blendingFactor::calcScheme
+(
+    const GeometricField<Type, fvPatchField, volMesh>& field,
+    const typename fv::convectionScheme<Type>& cs
+)
 {
-    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
-
-    if (!obr_.foundObject<fieldType>(fieldName_))
+    if (!isA<fv::gaussConvectionScheme<Type>>(cs))
     {
+        WarningInFunction
+            << "Scheme for field " << field.name() << " is not a "
+            << fv::gaussConvectionScheme<Type>::typeName
+            << " scheme. Not calculating " << resultName_ << endl;
+
         return;
     }
 
-    const fvMesh& mesh = refCast<const fvMesh>(obr_);
-
-    const fieldType& field = mesh.lookupObject<fieldType>(fieldName_);
-
-    const word divScheme("div(" + phiName_ + ',' + fieldName_ + ')');
-    ITstream& its = mesh.divScheme(divScheme);
-
-    const surfaceScalarField& phi =
-        mesh.lookupObject<surfaceScalarField>(phiName_);
-
-    tmp<fv::convectionScheme<Type>> cs =
-        fv::convectionScheme<Type>::New(mesh, phi, its);
-
     const fv::gaussConvectionScheme<Type>& gcs =
-        refCast<const fv::gaussConvectionScheme<Type>>(cs());
+        refCast<const fv::gaussConvectionScheme<Type>>(cs);
 
-    const surfaceInterpolationScheme<Type>& interpScheme =
-        gcs.interpScheme();
+
+    const surfaceInterpolationScheme<Type>& interpScheme = gcs.interpScheme();
 
     if (!isA<blendedSchemeBase<Type>>(interpScheme))
     {
-        FatalErrorInFunction
-            << interpScheme.typeName << " is not a blended scheme"
-            << exit(FatalError);
+        WarningInFunction
+            << interpScheme.type() << " is not a blended scheme"
+            << ". Not calculating " << resultName_ << endl;
+
+        return;
     }
 
     // Retrieve the face-based blending factor
@@ -123,4 +120,46 @@ void Foam::blendingFactor::calc()
 }
 
 
+template<class Type>
+void Foam::blendingFactor::calc()
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+
+    if (!obr_.foundObject<fieldType>(fieldName_))
+    {
+        return;
+    }
+
+    const fvMesh& mesh = refCast<const fvMesh>(obr_);
+
+    const fieldType& field = mesh.lookupObject<fieldType>(fieldName_);
+
+    const word divScheme("div(" + phiName_ + ',' + fieldName_ + ')');
+    ITstream& its = mesh.divScheme(divScheme);
+
+    const surfaceScalarField& phi =
+        mesh.lookupObject<surfaceScalarField>(phiName_);
+
+    tmp<fv::convectionScheme<Type>> tcs =
+        fv::convectionScheme<Type>::New(mesh, phi, its);
+
+    const fv::convectionScheme<Type>& cs = tcs();
+
+    if (isA<fv::boundedConvectionScheme<Type>>(cs))
+    {
+        const fv::boundedConvectionScheme<Type>& bcs =
+            refCast<const fv::boundedConvectionScheme<Type>>(cs);
+
+        calcScheme(field, bcs.scheme());
+    }
+    else
+    {
+        const fv::gaussConvectionScheme<Type>& gcs =
+            refCast<const fv::gaussConvectionScheme<Type>>(cs);
+
+        calcScheme(field, gcs);
+    }
+}
+
+
 // ************************************************************************* //

From 582ee21ef0d5ad60a07bc6e36437533ee4716c60 Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Thu, 11 Aug 2016 12:14:34 +0100
Subject: [PATCH 14/26] ENH: plane - updated construct from components to
 optionally normalise the normal vector

---
 .../meshes/primitiveShapes/plane/plane.C      | 20 ++++++++++++-------
 .../meshes/primitiveShapes/plane/plane.H      |  7 ++++++-
 2 files changed, 19 insertions(+), 8 deletions(-)

diff --git a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C
index a6d05b7f73..5f22c04431 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C
+++ b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C
@@ -129,23 +129,29 @@ Foam::plane::plane(const vector& normalVector)
 }
 
 
-Foam::plane::plane(const point& basePoint, const vector& normalVector)
+Foam::plane::plane
+(
+    const point& basePoint,
+    const vector& normalVector,
+    const bool normalise
+)
 :
     unitVector_(normalVector),
     basePoint_(basePoint)
 {
-    scalar magUnitVector(mag(unitVector_));
+    scalar magSqrUnitVector(magSqr(unitVector_));
 
-    if (magUnitVector > VSMALL)
-    {
-        unitVector_ /= magUnitVector;
-    }
-    else
+    if (magSqrUnitVector < VSMALL)
     {
         FatalErrorInFunction
             << "plane normal has zero length. basePoint:" << basePoint_
             << abort(FatalError);
     }
+
+    if (normalise)
+    {
+        unitVector_ /= Foam::sqrt(magSqrUnitVector);
+    }
 }
 
 
diff --git a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.H b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.H
index 44564fa92a..6645fb1b95 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.H
@@ -131,7 +131,12 @@ public:
         plane(const vector& normalVector);
 
         //- Construct from normal vector and point in plane
-        plane(const point& basePoint, const vector& normalVector);
+        plane
+        (
+            const point& basePoint,
+            const vector& normalVector,
+            const bool normalise = true
+        );
 
         //- Construct from three points in plane
         plane(const point& point1, const point& point2, const point& point3);

From 6afc9fb15d7e8a4ba919f6667778b6a2a2534158 Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Thu, 11 Aug 2016 12:14:43 +0100
Subject: [PATCH 15/26] ENH: tetOverlapVolume - improved robustness.  Fixes
 #207

---
 .../tetOverlapVolume/tetOverlapVolume.C       |   6 -
 .../tetOverlapVolume/tetOverlapVolume.H       |   6 -
 .../tetOverlapVolumeTemplates.C               | 105 ++++++++++++------
 3 files changed, 72 insertions(+), 45 deletions(-)

diff --git a/src/meshTools/tetOverlapVolume/tetOverlapVolume.C b/src/meshTools/tetOverlapVolume/tetOverlapVolume.C
index 209c2f7e37..cc54720dd3 100644
--- a/src/meshTools/tetOverlapVolume/tetOverlapVolume.C
+++ b/src/meshTools/tetOverlapVolume/tetOverlapVolume.C
@@ -39,12 +39,6 @@ namespace Foam
     defineTypeNameAndDebug(tetOverlapVolume, 0);
 }
 
-// When to consider a tet to be zero volume. We want to avoid doing clipping
-// against negative volume tets. Tet volume can be calculated incorrectly
-// due to truncation errors. The value below works for single and double
-// precision but could probably be tighter for double precision.
-Foam::scalar Foam::tetOverlapVolume::minTetVolume_ = SMALL*SMALL;
-
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
diff --git a/src/meshTools/tetOverlapVolume/tetOverlapVolume.H b/src/meshTools/tetOverlapVolume/tetOverlapVolume.H
index b027599218..a77e608c6e 100644
--- a/src/meshTools/tetOverlapVolume/tetOverlapVolume.H
+++ b/src/meshTools/tetOverlapVolume/tetOverlapVolume.H
@@ -54,12 +54,6 @@ class polyMesh;
 
 class tetOverlapVolume
 {
-    // Private data
-
-        //- Minimum tet volume to skip test
-        static scalar minTetVolume_;
-
-
     // Private classes
 
         //- tetPoints handling : sum resulting volumes
diff --git a/src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C b/src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C
index 3ff950311f..34eded1ff3 100644
--- a/src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C
+++ b/src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,51 +45,90 @@ void Foam::tetOverlapVolume::tetTetOverlap
     tetPointRef::storeOp cutInside(cutInsideTets, nCutInside);
     tetPointRef::dummyOp outside;
 
-    if (tetA.tet().mag() < minTetVolume_ || tetB.tet().mag() < minTetVolume_)
+    // Precompute the tet face areas and exit early if any face area is
+    // too small
+    static FixedList<vector, 4> tetAFaceAreas;
+    static FixedList<scalar, 4> tetAMag2FaceAreas;
+    tetPointRef tetATet = tetA.tet();
+    for (label facei = 0; facei < 4; ++facei)
     {
-        return;
+        tetAFaceAreas[facei] = -tetATet.tri(facei).normal();
+        tetAMag2FaceAreas[facei] = magSqr(tetAFaceAreas[facei]);
+        if (tetAMag2FaceAreas[facei] < ROOTVSMALL)
+        {
+            return;
+        }
     }
 
-    // face0
-    plane pl0(tetB[1], tetB[3], tetB[2]);
-    tetA.tet().sliceWithPlane(pl0, cutInside, outside);
-    if (nCutInside == 0)
+    static FixedList<vector, 4> tetBFaceAreas;
+    static FixedList<scalar, 4> tetBMag2FaceAreas;
+    tetPointRef tetBTet = tetB.tet();
+    for (label facei = 0; facei < 4; ++facei)
     {
-        return;
+        tetBFaceAreas[facei] = -tetBTet.tri(facei).normal();
+        tetBMag2FaceAreas[facei] = magSqr(tetBFaceAreas[facei]);
+        if (tetBMag2FaceAreas[facei] < ROOTVSMALL)
+        {
+            return;
+        }
     }
 
-    // face1
-    plane pl1(tetB[0], tetB[2], tetB[3]);
-    nInside = 0;
-    for (label i = 0; i < nCutInside; i++)
+
+    // Face 0
     {
-        const tetPointRef t = cutInsideTets[i].tet();
-        t.sliceWithPlane(pl1, inside, outside);
-    }
-    if (nInside == 0)
-    {
-        return;
+        vector n = tetBFaceAreas[0]/Foam::sqrt(tetBMag2FaceAreas[0]);
+        plane pl0(tetBTet.tri(0).a(), n, false);
+
+        tetA.tet().sliceWithPlane(pl0, cutInside, outside);
+        if (nCutInside == 0)
+        {
+            return;
+        }
     }
 
-    // face2
-    plane pl2(tetB[0], tetB[3], tetB[1]);
-    nCutInside = 0;
-    for (label i = 0; i < nInside; i++)
+    // Face 1
     {
-        const tetPointRef t = insideTets[i].tet();
-        t.sliceWithPlane(pl2, cutInside, outside);
-    }
-    if (nCutInside == 0)
-    {
-        return;
+        vector n = tetBFaceAreas[1]/Foam::sqrt(tetBMag2FaceAreas[1]);
+        plane pl1(tetBTet.tri(1).a(), n, false);
+
+        nInside = 0;
+        for (label i = 0; i < nCutInside; i++)
+        {
+            const tetPointRef t = cutInsideTets[i].tet();
+            t.sliceWithPlane(pl1, inside, outside);
+        }
+        if (nInside == 0)
+        {
+            return;
+        }
     }
 
-    // face3
-    plane pl3(tetB[0], tetB[1], tetB[2]);
-    for (label i = 0; i < nCutInside; i++)
+    // Face 2
     {
-        const tetPointRef t = cutInsideTets[i].tet();
-        t.sliceWithPlane(pl3, insideOp, outside);
+        vector n = tetBFaceAreas[2]/Foam::sqrt(tetBMag2FaceAreas[2]);
+        plane pl2(tetBTet.tri(2).a(), n, false);
+
+        nCutInside = 0;
+        for (label i = 0; i < nInside; i++)
+        {
+            const tetPointRef t = insideTets[i].tet();
+            t.sliceWithPlane(pl2, cutInside, outside);
+        }
+        if (nCutInside == 0)
+        {
+            return;
+        }
+    }
+
+    // Face 3
+    {
+        vector n = tetBFaceAreas[3]/Foam::sqrt(tetBMag2FaceAreas[3]);
+        plane pl3(tetBTet.tri(3).a(), n, false);
+        for (label i = 0; i < nCutInside; i++)
+        {
+            const tetPointRef t = cutInsideTets[i].tet();
+            t.sliceWithPlane(pl3, insideOp, outside);
+        }
     }
 }
 

From 9fb439bb46bdc5b488a3a8ccbc5b7209b6da4518 Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Thu, 18 Aug 2016 11:22:25 +0100
Subject: [PATCH 16/26] BUG: Corrected FFTW_ARCH_PATH - patch contributed by
 Bruno Santos.  Fixes #210

---
 etc/config.sh/FFTW | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/etc/config.sh/FFTW b/etc/config.sh/FFTW
index 7c0a48861d..f3cb4ea70e 100644
--- a/etc/config.sh/FFTW
+++ b/etc/config.sh/FFTW
@@ -64,7 +64,7 @@ then
     # it is either located within ThirdParty, or a central installation
     # outside of ThirdParty and must be added to the lib-path.
 
-    ending="${FFTW_ARCH_PATH_PATH##*-}"
+    ending="${FFTW_ARCH_PATH##*-}"
     if [ "$ending" != none -a "$ending" != system ]
     then
         _foamAddLib $FFTW_ARCH_PATH/lib$WM_COMPILER_LIB_ARCH

From 16b516bba3cf2590a482b493e17524fd6d30ddb7 Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Thu, 18 Aug 2016 15:24:13 +0100
Subject: [PATCH 17/26] ENH: noise functionality robustness improvements

---
 src/randomProcesses/noise/noiseFFT/noiseFFT.C   |  9 ++++++---
 .../noiseModels/surfaceNoise/surfaceNoise.C     | 17 +++++++++++++++--
 2 files changed, 21 insertions(+), 5 deletions(-)

diff --git a/src/randomProcesses/noise/noiseFFT/noiseFFT.C b/src/randomProcesses/noise/noiseFFT/noiseFFT.C
index 3a61c9b0dd..e82d25b90f 100644
--- a/src/randomProcesses/noise/noiseFFT/noiseFFT.C
+++ b/src/randomProcesses/noise/noiseFFT/noiseFFT.C
@@ -127,10 +127,13 @@ void Foam::noiseFFT::octaveBandInfo
 
     fBandIDs = bandIDs.toc();
 
-    // Remove the last centre frequency (beyond upper frequency limit)
-    fc.remove();
+    if (fc.size())
+    {
+        // Remove the last centre frequency (beyond upper frequency limit)
+        fc.remove();
 
-    fCentre.transfer(fc);
+        fCentre.transfer(fc);
+    }
 }
 
 
diff --git a/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.C b/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.C
index b8e312c0ce..92a3df57fd 100644
--- a/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.C
+++ b/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.C
@@ -494,8 +494,21 @@ void surfaceNoise::calculate()
         octave13FreqCentre
     );
 
-    List<scalarField> surfPSD13f(octave13BandIDs.size() - 1);
-    List<scalarField> surfPrms13f2(octave13BandIDs.size() - 1);
+    label bandSize = 0;
+    if (octave13BandIDs.empty())
+    {
+        WarningInFunction
+            << "Ocatve band calculation failed (zero sized). "
+            << "please check your input data"
+            << endl;
+    }
+    else
+    {
+        bandSize = octave13BandIDs.size() - 1;
+    }
+
+    List<scalarField> surfPSD13f(bandSize);
+    List<scalarField> surfPrms13f2(bandSize);
     forAll(surfPSD13f, freqI)
     {
         surfPSD13f[freqI].setSize(nLocalFace);

From 8e3dd0942556a21d8710b50ef4ffe3abe4301a31 Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Thu, 18 Aug 2016 15:29:56 +0100
Subject: [PATCH 18/26] BUG: ensightSurfaceReader - enabled reading of data
 from sub-directory and updated field mask (See #215)

---
 .../readers/ensight/ensightSurfaceReader.C    | 112 ++++++++++++------
 .../readers/ensight/ensightSurfaceReader.H    |  33 ++++--
 .../ensight/ensightSurfaceReaderTemplates.C   |  32 ++++-
 3 files changed, 128 insertions(+), 49 deletions(-)

diff --git a/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReader.C b/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReader.C
index 3e1b422bde..f90d05b85d 100644
--- a/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReader.C
+++ b/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReader.C
@@ -34,7 +34,7 @@ namespace Foam
 }
 
 
-void Foam::ensightSurfaceReader::skip(const label n, IFstream& is) const
+void Foam::ensightSurfaceReader::skip(const label n, Istream& is) const
 {
     label i = 0;
     token t;
@@ -59,6 +59,40 @@ void Foam::ensightSurfaceReader::skip(const label n, IFstream& is) const
 }
 
 
+void Foam::ensightSurfaceReader::readLine(IFstream& is, string& buffer) const
+{
+    buffer = "";
+    while (is.good() && buffer == "")
+    {
+        is.getLine(buffer);
+    }
+}
+
+
+void Foam::ensightSurfaceReader::debugSection
+(
+    const word& expected,
+    IFstream& is
+) const
+{
+    string actual = "";
+    readLine(is, actual);
+
+    if (expected != actual)
+    {
+        FatalIOErrorInFunction(is)
+            << "Expected section header '" << expected
+            << "' but read the word '" << actual << "'"
+            << exit(FatalIOError);
+    }
+
+    if (debug)
+    {
+        Info<< "Read section header: " << expected << endl;
+    }
+}
+
+
 void Foam::ensightSurfaceReader::readGeometryHeader(ensightReadFile& is) const
 {
     // Binary flag string if applicable
@@ -101,29 +135,6 @@ void Foam::ensightSurfaceReader::readGeometryHeader(ensightReadFile& is) const
 }
 
 
-void Foam::ensightSurfaceReader::debugSection
-(
-    const word& expected,
-    IFstream& is
-) const
-{
-    word actual(is);
-
-    if (expected != actual)
-    {
-        FatalIOErrorInFunction(is)
-            << "Expected section header '" << expected
-            << "' but read the word '" << actual << "'"
-            << exit(FatalIOError);
-    }
-
-    if (debug)
-    {
-        Info<< "Read section header: " << expected << endl;
-    }
-}
-
-
 void Foam::ensightSurfaceReader::readCase(IFstream& is)
 {
     if (debug)
@@ -141,31 +152,54 @@ void Foam::ensightSurfaceReader::readCase(IFstream& is)
     string buffer;
 
     // Read the file
+
     debugSection("FORMAT", is);
-    skip(3, is); // type: ensight gold
+    readLine(is, buffer);  // type: ensight gold
 
     debugSection("GEOMETRY", is);
-    readSkip(is, 2, meshFileName_);
+    readLine(is, buffer);
+    readFromLine(2, buffer, meshFileName_); // model: 1 xxx.0000.mesh
 
     debugSection("VARIABLE", is);
 
+    // Read the field description
     DynamicList<word> fieldNames(10);
-    DynamicList<word> fieldFileNames(10);
+    DynamicList<string> fieldFileNames(10);
     word fieldName;
-    word fieldFileName;
+    string fieldFileName;
     while (is.good())
     {
-        word primitiveType(is); // scalar, vector
+        readLine(is, buffer);
 
-        if (primitiveType == "TIME")
+        if (buffer == "TIME")
         {
             break;
         }
 
-        readSkip(is, 3, fieldName); // p, U etc
+        IStringStream iss(buffer);
+
+        // Read the field name, e.g. p U etc
+        readFromLine(4, iss, fieldName);
         fieldNames.append(fieldName);
 
-        is >> fieldFileName; // surfaceName.****.fieldName
+        // Field file name may contain /'s e.g.
+        //   surfaceName.****.fieldName
+        // This is not parser friendly - simply take remainder of buffer
+        label iPos = iss.stdStream().tellg();
+        fieldFileName = buffer(iPos, buffer.size() - iPos);
+        size_t p0  = fieldFileName.find_first_not_of(' ');
+        if (p0 == string::npos)
+        {
+            WarningInFunction
+                << "Error reading field file name. "
+                << "Current buffer: " << buffer
+                << endl;
+        }
+        else
+        {
+            size_t p1  = fieldFileName.find_last_not_of(' ');
+            fieldFileName = fieldFileName.substr(p0, p1 - p0 + 1);
+        }
         fieldFileNames.append(fieldFileName);
     }
     fieldNames_.transfer(fieldNames);
@@ -178,10 +212,14 @@ void Foam::ensightSurfaceReader::readCase(IFstream& is)
     }
 
     // Start reading time information
-    skip(3, is); // time set: 1
-    readSkip(is, 3, nTimeSteps_);
-    readSkip(is, 3, timeStartIndex_);
-    readSkip(is, 2, timeIncrement_);
+    readLine(is, buffer); // time set: 1
+
+    readLine(is, buffer);
+    readFromLine(3, buffer, nTimeSteps_);
+    readLine(is, buffer);
+    readFromLine(3, buffer, timeStartIndex_);
+    readLine(is, buffer);
+    readFromLine(2, buffer, timeIncrement_);
 
     if (debug)
     {
@@ -191,7 +229,7 @@ void Foam::ensightSurfaceReader::readCase(IFstream& is)
     }
 
     // Read the time values
-    skip(2, is);
+    readLine(is, buffer); // time values:
     timeValues_.setSize(nTimeSteps_);
     for (label i = 0; i < nTimeSteps_; i++)
     {
diff --git a/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReader.H b/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReader.H
index 8c99cfd466..0ba7e1c49b 100644
--- a/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReader.H
+++ b/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReader.H
@@ -69,7 +69,7 @@ protected:
         List<word> fieldNames_;
 
         //- Field file names
-        List<word> fieldFileNames_;
+        List<string> fieldFileNames_;
 
         //- Number of time steps
         label nTimeSteps_;
@@ -91,21 +91,38 @@ protected:
 
     // Protected Member Functions
 
+        //- Helper function to skip forward n steps in stream
+        void skip(const label n, Istream& is) const;
+
+        //- Helper function to read an ascii line from file
+        void readLine(IFstream& is, string& buffer) const;
+
         //- Read and check a section header
         void debugSection(const word& expected, IFstream& is) const;
 
-        //- Read the case file
-        void readCase(IFstream& is);
-
-        //- Helper function to skip forward n steps in stream
-        void skip(const label n, IFstream& is) const;
-
         //- Read (and throw away) geometry file header
         void readGeometryHeader(ensightReadFile& is) const;
 
+        //- Read the case file
+        void readCase(IFstream& is);
+
         //- Helper function to return Type after skipping n tokens
         template<class Type>
-        void readSkip(IFstream& is, const label nSkip, Type& value) const;
+        void readFromLine
+        (
+            const label nSkip,
+            IStringStream& is,
+            Type& value
+        ) const;
+
+        //- Helper function to return Type after skipping n tokens
+        template<class Type>
+        void readFromLine
+        (
+            const label nSkip,
+            const string& buffer,
+            Type& value
+        ) const;
 
         //- Helper function to return a field
         template<class Type>
diff --git a/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReaderTemplates.C b/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReaderTemplates.C
index 04a593575a..dace070320 100644
--- a/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReaderTemplates.C
+++ b/src/sampling/sampledSurface/readers/ensight/ensightSurfaceReaderTemplates.C
@@ -29,10 +29,10 @@ License
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 template<class Type>
-void Foam::ensightSurfaceReader::readSkip
+void Foam::ensightSurfaceReader::readFromLine
 (
-    IFstream& is,
     const label nSkip,
+    IStringStream& is,
     Type& value
 ) const
 {
@@ -42,6 +42,20 @@ void Foam::ensightSurfaceReader::readSkip
 }
 
 
+template<class Type>
+void Foam::ensightSurfaceReader::readFromLine
+(
+    const label nSkip,
+    const string& buffer,
+    Type& value
+) const
+{
+    IStringStream is(buffer);
+
+    readFromLine(nSkip, is, value);
+}
+
+
 template<class Type>
 Foam::tmp<Foam::Field<Type> > Foam::ensightSurfaceReader::readField
 (
@@ -60,9 +74,19 @@ Foam::tmp<Foam::Field<Type> > Foam::ensightSurfaceReader::readField
     fileName fieldFileName(fieldFileNames_[fieldIndex]);
 
     std::ostringstream oss;
-    oss << std::setfill('0') << std::setw(4) << fileIndex;
+    label nMask = 0;
+    for (size_t chari = 0; chari < fieldFileName.size(); chari++)
+    {
+        if (fieldFileName[chari] == '*')
+        {
+            nMask++;
+        }
+    }
+
+    const std::string maskStr(nMask, '*');
+    oss << std::setfill('0') << std::setw(nMask) << fileIndex;
     const word indexStr = oss.str();
-    fieldFileName.replace("****", indexStr);
+    fieldFileName.replace(maskStr, indexStr);
 
     
     ensightReadFile is(baseDir_/fieldFileName, streamFormat_);

From ac42dd2494944a757c39c4e3146298824032df30 Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Thu, 25 Aug 2016 13:32:45 +0100
Subject: [PATCH 19/26] BUG: mesh/filter case - corrected Allrun script.  Fixes
 #219

---
 tutorials/mesh/parallel/filter/Allrun | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/tutorials/mesh/parallel/filter/Allrun b/tutorials/mesh/parallel/filter/Allrun
index 8a0978f43f..4be640030b 100755
--- a/tutorials/mesh/parallel/filter/Allrun
+++ b/tutorials/mesh/parallel/filter/Allrun
@@ -4,6 +4,8 @@ cd ${0%/*} || exit 1    # Run from this directory
 # Source tutorial run functions
 . $WM_PROJECT_DIR/bin/tools/RunFunctions
 
+application=$(getApplication)
+
 # Create mesh
 runApplication blockMesh
 
@@ -16,7 +18,7 @@ runApplication topoSet
 # Create baffles and fields
 runApplication createBaffles -overwrite
 
-runApplication $(getApplication)
+runApplication $application
 
 #- RedistributePar to do decomposition
 runParallel redistributePar -decompose -cellDist

From bc5c6d83dfee70c62bf4afe34ce48c5b54fdc16d Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Thu, 25 Aug 2016 13:34:17 +0100
Subject: [PATCH 20/26] BUG: incompressible/icoFoam/cavityMappingTest - renamed
 Allrun-serial->Allrun.  Fixes #220

---
 .../icoFoam/cavityMappingTest/{Allrun-serial => Allrun}           | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename tutorials/incompressible/icoFoam/cavityMappingTest/{Allrun-serial => Allrun} (100%)

diff --git a/tutorials/incompressible/icoFoam/cavityMappingTest/Allrun-serial b/tutorials/incompressible/icoFoam/cavityMappingTest/Allrun
similarity index 100%
rename from tutorials/incompressible/icoFoam/cavityMappingTest/Allrun-serial
rename to tutorials/incompressible/icoFoam/cavityMappingTest/Allrun

From 16e4b25ac09636fe8de64b616e64c90fefe00395 Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Thu, 25 Aug 2016 13:36:43 +0100
Subject: [PATCH 21/26] BUG: Corrected Doxygen style sheet - patch contributed
 by Bruno Santos.  Fixes #217

---
 doc/Doxygen/css/doxyMod.css | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/doc/Doxygen/css/doxyMod.css b/doc/Doxygen/css/doxyMod.css
index 46ac2dfc5e..fbb968f594 100644
--- a/doc/Doxygen/css/doxyMod.css
+++ b/doc/Doxygen/css/doxyMod.css
@@ -118,6 +118,8 @@ tr.memlist
 .OFPlainTable tr td {
     height: 20px;
     padding-left: 5px;
+}
+
 div.line,
 span.comment,
 span.keyword,

From 6ffdce83aa460663f6249cf2927c8fb44efaad5d Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Thu, 25 Aug 2016 17:28:50 +0100
Subject: [PATCH 22/26] STYLE: Corrected spelling error

---
 src/thermophysicalModels/doc/thermophysicalModels.dox | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/thermophysicalModels/doc/thermophysicalModels.dox b/src/thermophysicalModels/doc/thermophysicalModels.dox
index a8f7f48ae9..9df1640a0b 100644
--- a/src/thermophysicalModels/doc/thermophysicalModels.dox
+++ b/src/thermophysicalModels/doc/thermophysicalModels.dox
@@ -23,7 +23,7 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-\page pageThermophsyicalModels Thermophsyical Models
+\page pageThermophsyicalModels Thermophysical Models
 
 \section secSchemes Overview
 The available thermophysical models are grouped into the following categories:

From 1589e8d431f563a61c3919f4a0660ea409d5049c Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Fri, 2 Sep 2016 08:14:10 +0100
Subject: [PATCH 23/26] BUG: buoyantPimpleFoam - div term should be present for
 both compressible and incompressible

---
 applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
index cda60d98b2..c778ef40a8 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
@@ -41,7 +41,6 @@
         p_rghDDtEqn =
         (
             fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
-          + fvc::div(phiHbyA)
         ==
             fvOptions(psi, p_rgh, rho.name())
         );
@@ -52,6 +51,7 @@
         fvScalarMatrix p_rghEqn
         (
             p_rghDDtEqn()
+          + fvc::div(phiHbyA)
           - fvm::laplacian(rhorAUf, p_rgh)
         );
 

From a4b62e6549c34eb6f809dda4b10cfcee5dc21fb1 Mon Sep 17 00:00:00 2001
From: sergio <sergio>
Date: Fri, 2 Sep 2016 11:51:59 -0700
Subject: [PATCH 24/26] Adding area weighted average to U mean calculation for
 uniformJump option in fan BC

---
 .../fields/fvPatchFields/derived/fan/fanFvPatchFields.C        | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchFields.C b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchFields.C
index d2580b1b82..fd0b61bdaa 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchFields.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchFields.C
@@ -46,7 +46,8 @@ void Foam::fanFvPatchField<Foam::scalar>::calcFanJump()
         scalarField Un(max(phip/patch().magSf(), scalar(0)));
         if (uniformJump_)
         {
-             Un = gAverage(Un);
+            scalar area = gSum(patch().magSf());
+            Un = gSum(Un*patch().magSf())/area;
         }
 
         if (phi.dimensions() == dimDensity*dimVelocity*dimArea)

From 292c522c71587f1837dcbb1bf9a88484fc35944d Mon Sep 17 00:00:00 2001
From: sergio <sergio>
Date: Fri, 2 Sep 2016 12:28:20 -0700
Subject: [PATCH 25/26] ENH: Adding iso-plane histogram to
 regionSizeDistribution

---
 .../regionSizeDistribution.C                  | 121 +++++++++++++++++-
 .../regionSizeDistribution.H                  |  46 ++++++-
 2 files changed, 164 insertions(+), 3 deletions(-)

diff --git a/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C
index e74bfb2f69..1d8855d665 100644
--- a/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C
+++ b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -335,7 +335,8 @@ Foam::regionSizeDistribution::regionSizeDistribution
     active_(true),
     alphaName_(dict.lookup("field")),
     patchNames_(dict.lookup("patches")),
-    log_(true)
+    log_(true),
+    isoPlanes_(dict.lookupOrDefault<bool>("isoPlanes", false))
 {
     // Check if the available mesh is an fvMesh, otherwise deactivate
     if (isA<fvMesh>(obr_))
@@ -388,6 +389,17 @@ void Foam::regionSizeDistribution::read(const dictionary& dict)
                 << "Transforming all vectorFields with coordinate system "
                 << coordSysPtr_().name() << endl;
         }
+
+        if (isoPlanes_)
+        {
+             dict.lookup("origin") >> origin_;
+             dict.lookup("direction") >> direction_;
+             dict.lookup("maxDiameter") >> maxDiameter_;
+             dict.lookup("nDownstreamBins") >> nDownstreamBins_;
+             dict.lookup("maxDownstream") >> maxDownstream_;
+
+             direction_ /= mag(direction_);
+        }
     }
 }
 
@@ -691,6 +703,8 @@ void Foam::regionSizeDistribution::write()
         // (patchRegions) and background regions from maps.
         // Note that we have to use mesh volume (allRegionVolume) and not
         // allRegionAlphaVolume since background might not have alpha in it.
+        // Deleting regions where the volume-alpha-weighted is lower than
+        // threshold
         forAllIter(Map<scalar>, allRegionVolume, vIter)
         {
             label regionI = vIter.key();
@@ -698,6 +712,7 @@ void Foam::regionSizeDistribution::write()
             (
                 patchRegions.found(regionI)
              || vIter() >= maxDropletVol
+             || (allRegionAlphaVolume[regionI]/vIter() < threshold_)
             )
             {
                 allRegionVolume.erase(vIter);
@@ -735,6 +750,108 @@ void Foam::regionSizeDistribution::write()
                 )
             );
 
+            vectorField centroids(sortedVols.size(), vector::zero);
+
+            // Check if downstream bins are calculed
+            if (isoPlanes_)
+            {
+                vectorField alphaDistance
+                (
+                    (alpha.internalField()*mesh.V())
+                  * (mesh.C().internalField() - origin_)()
+                );
+
+                Map<vector> allRegionAlphaDistance
+                (
+                    regionSum
+                    (
+                        regions,
+                        alphaDistance
+                    )
+                );
+
+                // 2. centroid
+                vectorField sortedMoment
+                (
+                    extractData
+                    (
+                        sortedRegions,
+                        allRegionAlphaDistance
+                    )
+                );
+
+                centroids = sortedMoment/sortedVols + origin_;
+
+                // Bin according to centroid
+                scalarField distToPlane((centroids - origin_) & direction_);
+
+                vectorField radialDistToOrigin
+                (
+                    (centroids - origin_) - (distToPlane*direction_)
+                );
+
+                const scalar deltaX = maxDownstream_/nDownstreamBins_;
+                labelList downstreamIndices(distToPlane.size(), -1);
+                forAll(distToPlane, i)
+                {
+                    if
+                    (
+                        (mag(radialDistToOrigin[i]) < maxDiameter_)
+                     && (distToPlane[i] < maxDownstream_)
+                    )
+                    {
+                        downstreamIndices[i] = distToPlane[i]/deltaX;
+                    }
+                }
+
+                scalarField binDownCount(nDownstreamBins_, 0.0);
+                forAll(distToPlane, i)
+                {
+                    if (downstreamIndices[i] != -1)
+                    {
+                        binDownCount[downstreamIndices[i]] += 1.0;
+                    }
+                }
+
+                // Write
+                if (Pstream::master())
+                {
+                    // Construct mids of bins for plotting
+                    pointField xBin(nDownstreamBins_);
+
+                    scalar x = 0.5*deltaX;
+                    forAll(xBin, i)
+                    {
+                        xBin[i] = point(x, 0, 0);
+                        x += deltaX;
+                    }
+
+                    const coordSet coords("distance", "x", xBin, mag(xBin));
+                    writeGraph(coords, "isoPlanes", binDownCount);
+                }
+
+                // Write to screen
+                if (log_)
+                {
+                    Info<< "    Iso-planes Bins:" << endl;
+                    Info<< "    " << token::TAB << "Bin"
+                        << token::TAB << "Min distance"
+                        << token::TAB << "Count:"
+                        << endl;
+
+                    scalar delta = 0.0;
+                    forAll(binDownCount, binI)
+                    {
+                        Info<< "    " << token::TAB << binI
+                            << token::TAB << delta
+                            << token::TAB << binDownCount[binI] << endl;
+                        delta += deltaX;
+                    }
+                    Info<< endl;
+
+                }
+            }
+
             // Calculate the diameters
             scalarField sortedDiameters(sortedVols.size());
             forAll(sortedDiameters, i)
diff --git a/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H
index 29c0bc5ac1..657077366b 100644
--- a/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H
+++ b/src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -35,6 +35,8 @@ Description
     based on where the field is below the threshold value.  These
     regions ("droplets") can now be analysed.
 
+
+
     Regions:
     - print the regions connected to a user-defined set of patches.
       (in spray calculation these form the liquid core)
@@ -53,6 +55,9 @@ Description
     - write graph of sum, average and deviation of droplet volume per bin
     - write graph of sum, average and deviation of user-defined fields.  For
       volVectorFields these are those of the 3 components and the magnitude.
+    - (optional) write graph of histogram of centroids on iso planes
+      downstream of the injector determined by origin, direction and maxDiameter
+      up to maxDownstream
 
     Example of function object specification:
     \verbatim
@@ -76,6 +81,24 @@ Description
             e3              (0 1 1);
             e1              (1 0 0);
         }
+
+        // Optional downstream iso-plane bins.
+        isoPlanes       true;
+
+            // Plane normal and point definition
+            direction       (1 0 1);
+            origin          (1e-4 0 5e-4);
+
+            // Maximum diameter of the cylinder formed by the origin point
+            // and direction
+            maxDiameter     3e-4;
+
+            // Maximum downstream distance
+            maxDownstream   6e-4;
+
+            // Number of iso-plane bins
+            nDownstreamBins 20;
+
     }
     \endverbatim
 
@@ -177,6 +200,27 @@ class regionSizeDistribution
         //- Optional coordinate system
         autoPtr<coordinateSystem> coordSysPtr_;
 
+        // Optional extra definition of bins on planes downstream to the origin
+        // point and maximum diameter
+
+            //- Switch to enable iso-planes sampling
+            bool isoPlanes_;
+
+            //- Optional origin point
+            vector origin_;
+
+            //- Optional plane normal direction
+            vector direction_;
+
+            //- Optional maximum diameter on plane
+            scalar maxDiameter_;
+
+            //- Optional number of bins for
+            scalar nDownstreamBins_;
+
+            //- Optional maximum downstream coordinate from origin
+            scalar maxDownstream_;
+
 
     // Private Member Functions
 

From 4dcbb32d242fdf27983659c1776fee203d288cc8 Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Wed, 7 Sep 2016 10:02:02 +0100
Subject: [PATCH 26/26] BUG: runTimePostProcessing - corrected surface
 visualisation caused by regression update/bug

---
 .../fieldVisualisationBase.C                  | 24 +++++++++++++++++--
 1 file changed, 22 insertions(+), 2 deletions(-)

diff --git a/src/postProcessing/functionObjects/graphics/runTimePostProcessing/fieldVisualisationBase.C b/src/postProcessing/functionObjects/graphics/runTimePostProcessing/fieldVisualisationBase.C
index 9d163da19d..86fd9b0e73 100644
--- a/src/postProcessing/functionObjects/graphics/runTimePostProcessing/fieldVisualisationBase.C
+++ b/src/postProcessing/functionObjects/graphics/runTimePostProcessing/fieldVisualisationBase.C
@@ -269,9 +269,29 @@ void Foam::fieldVisualisationBase::setField
             lut->SetVectorMode(vtkScalarsToColors::MAGNITUDE);
 
             // Configure the mapper
-            mapper->SelectColorArray(colourFieldName.c_str());
+            const char* fieldName = colourFieldName.c_str();
+            mapper->SelectColorArray(fieldName);
+
+            // Set to use either point or cell data
+            // Note: if both point and cell data exists, preferentially
+            //       choosing point data.  This is often the case when using
+            //       glyphs
+            if (pData->GetPointData()->HasArray(fieldName) == 1)
+            {
+                mapper->SetScalarModeToUsePointFieldData();
+            }
+            else if (pData->GetCellData()->HasArray(fieldName) == 1)
+            {
+                mapper->SetScalarModeToUseCellFieldData();
+            }
+            else
+            {
+                WarningInFunction
+                    << "Unable to determine cell or point data type "
+                    << "- assuming point data";
+                mapper->SetScalarModeToUsePointFieldData();
+            }
             mapper->SetScalarRange(range_.first(), range_.second());
-            mapper->SetScalarModeToDefault(); // try points, then cells
             mapper->SetColorModeToMapScalars();
             mapper->SetLookupTable(lut);
             mapper->ScalarVisibilityOn();